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series;  and (5)  concrete  examples  from  the  application  of  Volterra 
series  method  to  nonlinear  circuit  problems. 

In  presenting  the  important  aspects  of  the  Volterra  series  method 
,  we  deal  with  the  applicability  and  the  convergence  of  the  series. 
Recursive  relationships  for  the  upper  bound  on  the  error  incurred  when 
the  Volterra  series  is  truncated  are  derived  for  an  important  class  of 
nonlinear  system  problems.  It  is  shown  how  this  upper  bound  can  be 
directly  related  to  the  1^  norm  of  the  linear  kernel  function  and  the 
bound  on  the  input  of  the  linearized  system.  An  example  of  determining 
the  bound  on  the  input  function  in  terms  of  the  1^  norm  of  the  linear 
kernel  function  to  assure  the  convergence  of  the  series  solution  is 
also  given. 

The  characterization  of  multiple  polynomial  type  nonlinearity  cir¬ 
cuits  (systems)  using  the  Volterra  series  is  done  by  determining  the 
Volterra  Kernels,  or  their  transforms.  An  approach  based  on  applying 
multi-dimensional  transforms  directly  to  a  set  of  system  equations  for 
systematically  obtaining  the  transform  domain  description  of  the  kernel 
functions  is  developed.  The  case  of  multiple  input  circuits  is  also 
treated. 


In  studying  the  response  of  nonlinear  circuits,  relationships 
between  sinusoidal  steady  state  and  zero-state  transient  response  and 
the  transform  domain  description  of  the  kernel  functions  is  given.  We 
make  use  of  the  "association  of  variables"  in  the  transform  domain  to 
obtain  the  zero-state  transient  response  and  show  how  a  symbolic 
analysis  on  a  digital  computer  will  provide  both  types  of  solutions  - 
namely,  sinusoidal  steady-state  and  zero-state  solutions. 

^  Algorithms  for  adapting  the  Volterra  series  method  for  the  com¬ 
puter-aided  distortion  and  zero-state  transient  analysis  of  mildly 
nonlinear  circuits  are  developed  next.  A  symbolic  approach  is  used 
which  not  only  improves  the  computational  efficiency  but  also  provides 
a  link  to  the  earlier  ideas  on  convergence,  etc.  An  overview  of  the 
digital  computer  program  PRANC,  which  uses  the  Volterra  series  method, 
is  also  given. 
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Analysis." 

ill 


TABLE  OF  CONTENTS 


Page 

LIST  OF  TABLES . vi  i 

LIST  OF  FIGURES . vi !  I 

LIST  OF  ABBREVIATIONS . 

CHAPTER  1  -  INTRODUCTION .  1 

1-1.  Introduction  ............  .  ......  1 

1-2.  Nonlinear  System  Analysis  Methods .  2 

1-3.  Objectives  of  the  Investigation .  6 

1- 4.  Organization  of  the  Report .  9 

CHAPTER  2  -  VOLTERRA  FUNCTIONAL  SERIES .  12 

2- 1.  Introduction .  12 

2-2.  Volterra  Functional  Series  .  13 

2-3.  Derivation  of  the  Volterra  Series:  An  Example  .  18 

2- 3.  Convergence  properties  of  Volterra  Series .  24 

CHAPTER  3  -  MULTI PLE-NODE,  MULTIPLE-NONLINEARITY  CIRCUIT  ANALYSIS  .  47 

3- 1.  Introduction . 47 

3-2.  Multi-dimensional  Transforms  .  49 

3-3.  A  Nonlinear  Differential  Equation .  56 

3-4.  Multiple-node/  Multiple-nonlinearity  Circuit  Analysis.  ...  64 

3-5.  Multiple  Input  Circuit  Analysis .  78 

3- 6.  An  Example  . .  83 

CHAPTER  4  -  STEADY-STATE  AND  TRANSIENT  ANALYSIS  .  90 

4- 1.  Introduction .  90 

4-2.  Sinusoidal  Steady-State  Analysis  ....  .  90 

4- 3.  Transient  Analysis  Using  Volterra  Series  .  98 

CHAPTER  5  -  COMPUTER-AIDED  ANALYSIS  USING  VOLTERRA  SERIES  .  112 

5- 1.  Introduction . 112 

5-2.  Why  a  Symbolic  Approach.  .  . . 114 

5-3.  Symbolic  Analysis  Method . 117 

5-4.  Spectrum  and  Distortion  Analysis  Algorithm  .  124 

5-5.  Time-Domain  Analysis  Algorithm  .  128 

5-6.  Program  PRANC. . 132 


iv 


CHAPTER  6  -  NONLINEAR  CIRCUIT  ANALYSIS  EXAMPLES  .  138 

6-1.  Introduction . 138 

6-2.  Spectrum  and  Distortion  Analysis  Example  .  139 

6-3.  Time-domain  Analysis  Example  .  147 

6- 4.  PRANC  vs.  NCAP . 149 

CHAPTER  7  -  CONCLUDING  REMARKS . 151 

7- 1.  Summary.  . . 151 

7-2.  Further  Research  .  .  ............  152 

REFERENCES . 154 

APPENDIX . 160 


J 


LIST  OF  TABLES 


TABLE  Page 

3- 1.  Nonlinear  Current  Sources  in  multiple-node, 

multiple-nonlinearity  circuit  analysis  .  73,74 

4- 1.  List  of  Associated  Transforms .  102 

6-1.  Transistor  Parameters  for  the  Circuit  of  Fig.  6-1 .  139 

6-2.  Coefficients  of  the  polynomial  type  nonlinearities  ....  143 

6-3.  Summary  of  Trunction  Error  .  149 


LIST  OF  FIGURES 


Figure  Page 

2-1.  A  Time-invariant  Nonlinear  System.  ..  .  13 

2-2.  Volterra  Series  Equivalent  System  Representation  .  17 

2-3.  A  System  with  "non-PS"  kernal .  35 

2- 4.  Roots  of  the  Nonlinear  Algebraic  Equation .  45 

3- 1.  Step  in  Nonlinear  Circuit  Analysis  Using  Volterra  Series  .  68,69 

3-2.  Determination  of  Volterra  Transfer  Functions  .  77 

3-3.  Multiple  Input  Nonlinear  Circuit  Analysis . 80,81 

3-4.  A  Nonlinear  Circuit  Example .  84 

5-1 .• Determination  of  Z(s)  for  the  p-port  network  using  state 

equations .  120 

5-2.  Algorithm  for  Inverting  Y<s)  symbolically .  123 

5-3.  Algorithm  for  Spectrum  and  Distortion  Analysis  .  126 

5- 4.  Algorithm  for  Computing  time-domain  response  using 

Volterra  series .  131 

6- 1.  A  Transistor  Amplifier  Circuit  .  140 

6-2.  Transistor  equivalent  model .  142 

6-3.  First-order  Transfer  function  magnitude  characteristics.  .  144 

6-4.  Second-order  Trap  -fer  function  magnitude  characteristics  .  144 

6-5.  Harmonic  Distortion  Variation  with  amplitude  .  146 

6-6.  A  Single  Nonlinearity  Circuit  Example .  147 


LIST  OF  ABBREVIATIONS 


FM  -  Frequency  Modulation 

PS  -  Product  Structure 

PFE  -  Partial  Fraction  Expansion 


CHAPTER  1 


INTRODUCTION 


1-1 .  Introduction 

The  superposition  property  shared  by  all  linear  systems  makes  it 
possible  to  develop  unified  methods  of  analysis  and  synthesis  of  these 
systems.  By  contrast,  nonlinear  systems  do  not  share  this  property  and 
as  such  no  single  method  exists  at  this  time  which  can  be  applied  to  the 
analysis  and  synthesis  of  all  such  systems. 

The  area  of  nonlinear  systems  analysis,  however,  remains  an  impor¬ 
tant  area  of  research.  Most  devices  which  are  linear  in  a  certain  re¬ 
gion  of  operation  may  reveal  a  nonlinear  behavior  when  operated  in  a 
different  region.  Devices  such  as  diodes,  transistors,  tunnel  diodes, 
etc.,  which  are  ubiquitous  to  all  electronic  circuits,  are  inherently 
nonlinear  (except  in  a  very  small  region  of  operation).  The  analysis  of 
circuits  containing  such  devices  is  an  important  problem. 

In  the  analysis  of  nonlinear  systems,  two  alternatives  are  avail¬ 
able:  The  first  possibility  is  to  take  a  problem-dependent  analysis  ap¬ 
proach;  that  is,  take  a  specific  nonlinear  system  and  perform  a  detailed 
analysis  on  it,  and,  in  the  process,  develop  some  efficient  method  of 
analysis.  In  general,  this  "new"  method  may  not  be  useful  for  other 
nonlinear  system  problems.  The  second  alternative  is  to  develop  and  ap¬ 
ply  an  analysis  method  which  is  applicable  to  a  broad  class  of  nonlinear 
systems,  and  Is  not  oriented  towards  a  specific  system. 
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In  the  present  day,  when  computer-aided  analysis  pervades  many 
areas  of  science  and  engineering,  the  computer  has  become  an  important 
arbiter  in  the  success  and  popularity  of  an  analysis  method.  Even  a 
fairly  general  analysis  method  may  remain  in  oblivion  if  it  is  found  - 
or  thought  -  unadaptable  to  computer-aided  analysis. 

The  Volterra  functional  series  HI ,21  for  nonlinear  system  analysis 
may  be  regarded  as  one  such  method:  it  is  applicable  to  a  large  class 
of  nonlinear  system  problems,  but  is  not  widely  used  in  engineering  cal¬ 
culations.  Although  much  theoretical  work  £3-6,12-14,15-241  has  been 
done  on  this  method  over  the  past  two  decades,  the  applications  of  this 
method  have  been  limited  in  their  scope.  Most  authors  of  these  publica¬ 
tions  have  ignored  the  computational  aspect  of  this  method,  and  those 
who  have  considered  it  have  deemed  this  method  as  cumbersome  or  unadapt¬ 
able  to  computer-aided  analysis.  In  this  report  we  will  use  the  Volter¬ 
ra  series  as  a  basis  for  computer-aided  circuit  analysis  of  nonlinear 
circuits,  and,  in  the  process,  introduce  some  interesting  aspects  of  the 
method. 

An  overview  of  some  nonlinear  system  analysis  methods  is  given  io 
section  1-2,  with  special  emphasis  on  the  Volterra  series  approach. 
Section  1-3  presents  the  objectives  of  this  investigation.  Finally,  a 
glimpse  of  the  results  of  this  research  effort,  along  with  the  organiza¬ 
tion  of  this  report,  is  presented  in  section  1-4. 

1-2.  Nonlinear  System  Analysis  Methods. 

It  is  indeed  rare  to  find  a  nonlinear  system  analysis  problem  which 
has  a  closed-form  solution.  One  must  then  resort  to  using  approximation 
methods  for  gaining  insight  into  the  system  behavior.  The  broad  ca- 
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tegory  of  approximation  methods  can  be  roughly  broken  into  two  groups: 
1)  analytical  methods,  which  often  yield  qualitative  and  quantitative 
information  about  the  system,  and  2)  numerical  methods,  which  give  quan¬ 
titive  information  about  the  system.  The  well-known  methods  belonging 
to  the  latter  group  are:  1)  polynomial  approximation  to  the  solution 
function,  leading  to  the  Numerical  Integration  methods;  and  2)  Taylor 
series  approximation  of  the  solution,  leading  to  the  Runge-Kutte  Type 
methods.  Our  survey  on  the  approximation  methods  will  deal  with  the 
analytical  methods  only.  We  consider  some  of  them  in  the  following 
sub-sections. 

J[-2.1_  Iteration  Method.  A  method  for  solving  nonlinear  differential  or 
integral  equations  is  based  on  the  process  of  successive  iteration, 
called  the  iteration  method.  This  process  can  be  performed  in  a  number 
of  ways,  but  the  basic  procedure  is  the  same.  The  nonlinear  equation  is 
first  solved  by  neglecting  certain  terms  -  generally  the  nonlinear  terms 
-  in  the  equation.  The  resulting  solution  is  then  re-substituted  in  the 
system  equation  without  neglecting  any  t?-rms.  This  process  is  repeated 
and,  under  well  defined  conditions  C15, 18,353,  each  resulting  solution 
is  a  better  approximation  to  the  actual  solution. 

When  nonlinear  differential/integral  equations  with  polynomial  type 
nonlinear  terms  are  solved  via  the  iteration  process,  the  implementation 
involves  the  actual  squaring,  cubing,  etc.,  operations.  When  sinusoidal 
steady  state  response  is  obtained  using  this  method,  the  approximation 
to  the  entire  output  spectrum  is  obtained  in  one  step,  in  contrast  to 
some  other  methods,  which  give  new  distortion  terms  with  each  step. 

1-2.2  Perturbation  Method.  Like  the  iteration  method,  the  perturbation 
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method  is  applied  in  a  wide  variety  of  ways.  This  method  is  generally 
applied  to  nonlinear  differential  equations,  in  which  a  small  parameter 
is  associated  with  the  nonlinear  terms,  by  introducing  a  dummy  variable 
to  help  dissolve  the  nonlinear  differential  equation  into  a  sequence 
linear  differential  equations  with  nonlinear  terms  which  can  be  solved 
in  a  bootstrapping  manner.  The  linearized  system  is  first  solved.  The 
resulting  solution  is  used  to  solve  the  next  system  equations  involving 
the  quadratic  terms.  This  is  followed  by  solving  the  linear  differen¬ 
tial  equations  with  cubic  terms,  and  so  on.  The  solutions  obtained  at 
each  step  are  added  to  yield  the  approximation  to  the  actual  solution, 
thus  giving  a  series  type  approximation. 

The  description  of  the  procedure  above  is  the  basic  perturbation 
method.  However,  in  some  nonlinear  problems  this  approach,  when  applied 
as  above,  leads  to  serious  convergence  problems.  The  method  is  then 
modified  to  get  rid  of  the  so-called  "secular  terms,"  which  grow  up  in¬ 
definitely  as  t  +  «  C17,183. 

1_-2.3.  Volterra  Series  Method.  The  Volterra  series  method  is  a  type  of 
functional  series  which  relates  the  system  input,  x(t),  to  the  system 
output,  y(t),  as: 

•  e  e  " 

y<t)  =  2  hn(T1,...,Tn)  n  x<t-T.)  dT.  (1-1) 
n=1  n-fold  n  1  n  i=1  11 

where  hn  is  the  n-th  order  kernel  function,  whose  n-dimensional 
transform  is  called  the  n-th  order  transfer  function.  Clearly,  the  sys¬ 
tem  characterization,  using  Volterha  series,  is  done  by  determining 
these  kernel  functions,  or  their  transforms.  The  solution  obtained  is 
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of  a  series  type. 

Wiener  C23  first  applied  this  method  to  calculate  the  response  of  a 
nonlinear  system  with  memory  to  a  Gaussian  noise  input.  Much  of 
theoretical  work  on  this  method  was  done  in  the  late  50's  and  early 
60's.  Brilliant  C43  studied  the  convergence  of  the  series  for  bounded 
inputs.  George  C53  devised  a  "system  of  algebra"  for  combining  Volterra 
systems  and  developed  the  "association  of  variables"  technique  for 
directly  going  from  a  multi-dimensional  transform  description  to  the 
one-dimensional  transform  description.  Zames  C33  studied  Volterra  sys¬ 
tems  when  placed  in  a  feedback  loop.  Bedrosian  and  Rice  C123  and  Rudko 
and  Wiener  C373  give  a  set  of  formulae  for  Volterra  systems  driven  by 
random  inputs.  This  set  of  formulae  have  been  applied  to  Gaussian  ran¬ 
dom  inputs  only. 

The  approach  has  also  been  applied  to  specific  problems.  Van  Trees 
C93  characterized  the  phase-locked  loop  with  a  nonlinearity.  Bedrosian 
and  Rice  C123  and  Baranyi  C383  used  the  method  to  calculate  and  compen¬ 
sate  the  distortion  incurred  in  an  FM  signal,  respectively.  Narayanan 
C73  and  Gopal,  et.  al.  C363  calculated  the  intermodulation  distortion  in 
transistor  amplifiers.  The  modeling  of  communication  receivers  using 
Volterra  series  was  done  in  C103.  Ewen  Cl 1 3  studied  the  identification 
of  Volterra  systems  and  Naditch  C133  applied  the  approach  for  high  fre¬ 
quency  calculations.  The  relationship  between  Volterra  series  and  Pi¬ 
card  iteration,  for  a  class  of  nonlinear  system  problems,  was  developed 
by  Leon  and  Shaefer  Cl 93. 
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1 -3 .  Objectives  of  the  Investigation. 

In  the  analysis  of  nonlinear  systems,  two  main  classes  of  solutions 
are  generally  sought:  1)  transient,  and  2)  steady  state.  The  basic 
goal  of  this  investigation  is  to  obtain  these  solutions  for  nonlinear 
circuits  via  the  Volterra  series  method.  Before  proceeding  with  the 
main  subject  of  this  sub-section,  we  briefly  provide  the  motivation  for 
using  Volterra  series. 

The  most  commonly  used  present-day  approach  for  analyzing  nonlinear 
systems  is  numerical  integration.  The  nonlinear  differential  equations 
are  integrated  from  some  initial  time,  tQ,  to  some  final  time,  t^.  When 
the  sinusoidal  steady  state  response  is  sought,  the  value  of  t^  chosen 
is  usually  Large  to  insure  that  all  transients  have  decayed.  A  fast 
Fourier  transform  then  yields  the  frequency  components  of  the  output 
response.  A  more  efficient  method  for  obtaining  the  sinusoidal  steady- 
state  response  is  to  pose  the  analysis  problem  as  a  two-point  boundary 
value  problem  C31,32]  and  then  apply  an  iterative  two  step  algorithm. 
This  approach,  however,  allows  only  single  frequency  inputs. 

The  problems  involved  in  the  numerical  integration  method  are  well 
known  C293.  These  problems  notwithstanding,  there  are  other  inefficien¬ 
cies.  When  one  is  solely  interested  in  the  sinusoidal  steady  state 
response,  the  computation  expended  in  reaching  close  to  t^  is  a  waste. 
This  inefficiency  grows  as  poles  of  the  linearized  system  get  close  to 
the  j  -axis,  as  is  often  the  case  in  many  quasi-linear  communication 
ci rcuits. 

Other  methods  such  as  the  harmonic  balance  or  the  describing  func¬ 
tion  method  are  seldom  used,  simply  because  the  assumption  behind  these 
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methods  render  them  undependable.  The  Picard  iteration  method  can  also 
be  used  for  nonlinear  system  analysis.  This  method  also  has  limitations 
when  used  for  computer-aided  analysis,  particularly  when  multi-tone  in¬ 
puts  are  present.  The  basic  perturbation  method  and  Volterra  series  are 
alike  when  applied  to  nonlinear  circuit  analysis  problems. 

We  now  return  to  the  main  subject  of  this  sub-section,  namely,  the 
goals  of  this  investigation. 

As  mentioned  in  the  previous  sub-section,  the  nonlinear  system  is 
completely  characterized  by  the  Volterra  kernel  functions,  or  their 
multi-dimensional  transforms,  referred  to  as  Volterra  transfer  func¬ 
tions.  In  analyzing  a  nonlinear  circuit,  the  first  problem  therefore  is 
to  characterize  it. 

Once  the  kernel  functions  or  the  transfer  functions  are  known,  we 
can  determine  the  response  of  the  circuit.  However,  since  the  Volterra 
series  is  an  infinite  series,  the  error  incurred  as  a  result  of  trunca¬ 
tion  at  a  finite  term  must  be  investigated.  This  leads  to  the  problem 
of  convergence.  We  look  into  how  the  series  converges  for  a  class  of 
lumped  nonlinear  circuits. 

Next  we  investigate  the  implementation  of  the  Volterra  series 
method  on  a  digital  computer.  In  doing  so,  we  first  develop  the  basic 
algorithm  for  adapting  this  method  for  computer-aided  analysis  and  then 
actually  implement  it  for  obtaining  the  steady-state  response  of  cir¬ 
cuits  with  multiple  nonlinearity,  multiple  multi-tone  input  sources. 
The  determination  of  the  transient  response  from  the  Volterra  transfer 
function  for  Laplace  transformable  inputs  is  also  investigated. 
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In  summary,  the  goals  of  this  investigation  can  be  categorically 
stated  as: 

1)  Devise  a  systematic  way  of  deriving  the  Volterra  kernel  func¬ 
tions,  or  their  transform  domain  description. 

2)  Determine  a  method  of  estimating  the  error  incurred  as  a  result 
of  truncating  the  series  solution  at  a  finite  term. 

3)  Obtain  an  algorithm  for  adapting  the  Volterra  series  method  for 
analyzing  multiple  nonlinearity  circuits. 

4)  Implement  the  algorithm  of  (3)  as  a  computer  program. 

Before  concluding  this  subsection,  we  make  the  following  comment: 
Many  papers  and  theses  on  the  theoretical  aspects  of  Volterra  series 
have  appeared  in  the  literature,  most  of  which  have  been  included  in  the 
References.  After  studying  these  references,  one  is  still  left  with  the 
most  practical  questions:  How  do  we  apply  this  method  to  the  real  world 
nonlinear  circuit  problems?  How  do  we  use  this  approach  for  analyzing 
and  computing  the  behavior  of  a  nonlinear  system  problem  at  hand? 

The  motivation  behind  this  research  was  to  seek  answers  to  the 
above  fundamental  questions  and  also  remove  the  misconception  that  the 
Volterra  series  method  is  cumbersome  for  computer  analysis.  We  feel 
that  this  method  could  be  employed  in  a  wider  variety  of  engineering 
calculations  than  is  presently  the  case.  The  contents  of  this  report 
are  therefore  centered  around  the  practical  aspects  concerning  Volterra 
series. 
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In  addition  to  this  introductory  chapter,  this  report  comprises  six 
chapters. 

Chapter  2,  entitled  "Volterra  Series",  discusses  the  analysis 
method  (1-2.3)  which  forms  the  basis  of  this  investigation.  Commonly 
used  terminology  and  the  salient  features  associated  with  the  series  are 
presented.  The  applicability  along  with  the  convergence  of  the  Volterra 
series  is  discussed.  In  treating  the  topic  of  convergence,  we  concern 
ourselves  not  with  abstract  nonlinear  equations,  but  with  a  concrete 
class  of  time-invariant  nonlinear  system  equations  with  polynomial  type 
nonlinearities.  Recursive  relationships,  which  provide  an  estimate  of 
the  error  due  to  truncating  the  series  (1-1),  for  this  class  of  non¬ 
linear  differential  equations  are  developed.  Knowing  the  ^  norm  of  the 
linear  kernel  function,  one  can  estimate  the  number  of  kernel  functions 
that  needs  be  derived  to  meet  a  prescribed  error  criterion.  When  the 
solution  is  sought  during  a  finite  time  interval,  the  conditions  for 
which  the  Volterra  series  converges  are  given  in  [21,223.  However,  for 
the  interval  convergence  of  the  series  has  not  been  proven.  We 
present  an  example,  based  on  C233  and  the  new  recursive  relationships, 
to  show  how  the  bound  on  the  input,  in  terms  of  the  norm  of  the 
linear  kernel  function,  for  which  the  series  in  this  interval  converges 
can  be  derived. 

In  Chapter  3,  we  deal  with  the  topic  of  multiple  nonlinearity  cir¬ 
cuit  analysis.  In  dealing  with  this  topic  we  develop  a  systematic 
method  of  obtaining  the  transform  domain  description  of  the  Volterra 
kernel  functions  for  a  large  class  of  nonlinear  systems.  Our  approach 


in  contrast 


relies  on  the  application  of  multi-dimensional  transforms, 
to  the  "harmonic  input"  method  used  by  other  authors  l!7,12-14].  This 
new  approach,  besides  being  more  intuitive,  avoids  the  morass  of  algebra 
involved  in  the  other  approaches.  The  determination  of  Volterra 
transfer  functions  for  multiple  nonlinearity  circuits  with  single  inputs 
is  described.  The  extension  to  multiple  inputs  is  also  discussed  brief¬ 
ly. 

Chapter  4  deals  with  the  determination  of  the  sinusoidal  steady- 
state  and  zero-state  transient  responses  from  the  Volterra  transfer 
functions.  The  case  of  multi-tone  inputs  is  considered,  and  the  deter¬ 
mination  of  the  responses  for  calculating  the  distortion  indices  at  the 
various  and  many  frequencies  is  discussed.  The  fact  that  convolution  in 
the  time-domain  becomes  multiplication  in  the  transform  domain,  along 
with  the  "association  of  variables"  technique,  is  exploited  to  show  how 
the  zero-state  transient  response  can  be  determined  for  inputs  which  are 
Laplace  transformable  and  factorable.  Nonlinear  lumped  systems  with 
such  inputs  obviate  the  need  for  any  numerical  integration.  With  the 
complete  pole-residue  information  of  the  linearized  system,  it  is  shown 
that  one  can  determine  the  time-domain  response  for  such  inputs. 

The  ideas  gathered  in  Chapters  2,  3,  and  4  are  used  to  develop  al¬ 
gorithms  for  most  "efficiently"  implementing  the  Volterra  series  method 
for  nonlinear  circuit  analysis  on  a  digital  computer.  This  is  the  topic 
of  Chapter  5.  The  fundamental  idea  used  in  this  development  is  the  use 
of  semi-symbolic  C29]  analysis  of  the  linearized  circuit.  Through  this 
approach,  many  advantages  are  accrued:  1)  repeated  inversion  of  a 
"large"  matrix  is  avoided;  2)  the  pole-residue  information  obtained  can 
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be  used  for  zero-state  transient  analysis;  3)  the  pole-residue  informa¬ 
tion  is  directly  related  to  the  l^  norm  of  the  linear  kernel  function, 
which,  in  turn,  is  related  to  the  truncation  error  (Chapter  2);  4)  the 

functional  description  of  the  higher  order  transfer  functions  can  be  ob¬ 
tained  from  the  functional  description  of  the  linear  functions  (Chapter 
3).  Thus,  Chapter  5,  in  a  sense,  lends  unity  to  this  whole  investiga¬ 
tion. 

Chapter  6  gives  an  example  from  the  use  of  PRANC.  Various  ramifi¬ 
cations  from  the  results  obtained  are  discussed. 

Chapter  7  is  reserved  for  concluding  remarks  and  some  topics  of 
further  research. 


II 


CHAPTER  2 


VOLTERRA  FUNCTIONAL  SERIES 


2-1 .  Introduction 

In  many  of  the  problems  encountered  in  nonlinear  systems  analysis, 
it  is  difficult  -  if  not  impossible  -  to  find  explicit  closed-form  ex¬ 
pressions  for  the  solution.  Approximation  for  the  solution  must  then  be 
sought.  Several  methods  [16-18]  which  yield  approximate  solution  are 
available,  and  were  discussed  in  the  previous  section. 

The  intention  behind  using  the  Volterra  functional  series  is  funda¬ 
mentally  to  evaluate  solutions  which  cannot  be  obtained  in  closed  form, 
and  therefore  the  solution  obtained  from  this  approach  will  involve  an 
approximation.  Whenever  approximation  methods  are  used,  the  question  of 
convergence  is  an  important  consideration.  In  this  section  we  present 
the  general  convergence  properties  of  the  Volterra  series  method.  In 
dealing  with  this  topic  we  rely  on  the  results  of  [22],  These  results 
are  extended  to  an  important  class  of  nonlinear  systems  by  establishing 
the  bounds  on  the  output  in  terms  of  the  bound  on  the  input  and  the  i ^ 
norm  of  the  impulse  response  of  the  linearized  system.  We  believe  that 
these  results  are  of  more  practical  value  to  a  nonlinear  system  analyst 
and  designer. 

Other  theoretical  considerations,  such  as  the  applicability  of  the 
method  and  input-output  uniqueness  are  also  briefly  discussed.  This 
discussion  relies  heavily  on  the  works  of  George  C53,  Brilliant  C4D. 

In  section  2-2  we  present  the  Volterra  functional  series  represen¬ 
tation  of  nonlinear  systems  along  with  the  applicability  of  the  method. 
In  discussing  the  applicability  we  present  a  sufficient  condition  C4] 


12 


for  a  nonlinear  system  to  be  representable  by  the  functional  series 
method. 

Section  2-3  is  used  to  illustrate  the  application  of  the  Volterra 
series  method  for  obtaining  the  approximate  solution  of  a  general  non¬ 
linear  time-invariant  differential  equation. 

Section  2-4  deals  with  the  convergence  properties  of  the  Volterra 
functional  series.  Both  the  single  nonlinearity  case  and  the  multiple 
nonlinearity  case  are  discussed  in  this  sub-section.  A  brief  exposition 
on  the  stability  of  solutions  is  also  presented.  The  stability  con¬ 
sidered  here  is  that  the  response,  y(t),  of  the  system  remains  bounded; 
that  is,  y(t)  <  <»,  for  all  t  _>  tg. 

2-2.  Volterra  Functional  Series 

Consider  the  time-invariant  nonlinear  system  of  Figure  2-1,  with  an 
input,  x(t),  and  an  output,  y(t>.  The  nonlinear  system  acts  as  an 
operator  that  maps  the  function  x(t)  into  a  function  y(t>;  that  is 

y (t)  =  SCx(t)3  (2-1) 


x(t) 


Fig.  2-1.  A  time-invariant  nonlinear  system 


13 


The  basic  idea  behind  the  VoLterra  functional  series  is  to  approximate  S 
by  a  series  of  functionals.  This  idea  is  analogous  to  the  case  of  func¬ 
tions,  where  any  function,  under  some  well-defined  conditions,  can  be 
approximated  by  a  series  of  functions,  say,  for  example,  polynomials. 
Furthermore,  for  polynomials,  the  Weierstrass  theorem  shows  that  any 
continuous  function  can  be  approximated  arbitrarily  closely  in  the  limit 
in  a  region  by  a  sequence  of  polynomials  C333.  Similarly,  it  has  been 
shown*  that  for  a  continuous  functional  SCx(t;  tg  <  t  <  t^)3,  there  ex¬ 
ists  a  sequence  of  functionals  which  approximate  S  arbitrarily  closely 
in  the- limit.  While  the  existence  of  these  sequence  of  functionals  is 
proved,  their  method  of  determination  is  not  given.**  The  Volterra  func¬ 
tional  series  is  therefore  a  type  of  functional  series  used  for  approxi¬ 
mating  S.  It  relates  the  input,  x(t),  to  the  output,  y(t>,  as  follows: 

y(t)  =  Z  //•••/  h  )  n  xCt-x  )dx  (2-2) 
n=1  n-fold  n  1  n  p=1  p  p 

where  h  (t.,...,t  )  is  known  as  the  time-domain  Volterra  kernel  of  order 
n  l  n 

n.  It  should  be  clear  that  the  nonlinear  system  characterization  is 
done  through  the  determination  of  these  kernel  functions;  for,  once  we 
know  them,  the  output  function  can  be  determined  from  eqn.  (2-2)  for  a 
known  input  function.  For  causal  kernel  functions,  the  limits  of  in¬ 
tegration  are  between  0  and  »,  i.e.  tg  is  taken  to  be  zero. 

*Van  Trees  C83  claims  that  this  was  done  by  Frechet.  For  original 
reference  see  C83. 

**This  again  is  analogous  to  the  case  of  functions:  the  existence  of 
the  sequence  of  approximating  functions  is  proved  by  Weierstrass, 
but  a  method  of  finding  these  functions  is  not  given.  We  then 
resort  to  using  orthogonal  polynomials,  etc.  as  basis  functions. 


Two  well-known  cases  are  readily  derived  from  eqn.  2-2:  First,  when 


the  set  of  kernels  th^h^,...}  are  zero,  eqn.  (2-2)  reduces  to  the 

well-known  convolution  integral  of  a  linear  system;  and,  two,  when  the 

n 

kernel  h  is  given  by  a  n  S(T  K  then  eqn.  (2-2)  represents  the  out- 

n  n  p=1  p 

put  of  a  nonlinear,  memoryless  system. 

Before  proceeding  with  the  discussion  of  the  applicability  of  Vol- 
terra  functional  series,  we  point  out  some  important  features  and  defin¬ 
itions  associated  with  the  series: 

(a)  The  term  order  with  respect  to  an  output  component  is  often 
used.  Order  in  this  context  is  defined  as  the  number  of  times  the  input 
is  multiplied  with  itself  in  the  convolution  integrals  appearing  in  eqn. 
2-2.  Due  to  this  definition  of  order,  it  is  often  more  convenient  to 
write  the  series  as: 


yCt)  =  £y  (t) 
n=1  n 


(2-3) 


where 


/•••/ 


n 


>  n  x(t-T  >dt 
P-1  p 


p 


(2-4) 


Equation  (2-3)  then  represents  the  output,  y(t),  as  a  summation  of  the 
outputs  of  order  1,2,3,.... 

(b)  As  shown  in  eqn.  (2-4),  the  nth  order  output  functional  is 

equal  to  the  convolution  of  the  kernel  It  with  the  n-fold  product  of  the 

n 

system  input  values.  Thus,  if  the  input  is  an  impulse  function,  then  hR 
can  be  viewed  as  the  nonlinear  impulse  response  function  of  order  n. 
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whose  n  dimensional  transform  can  then  be  called  the  nonlinear  transfer 
function  of  order  n  [10,143. 

(c)  The  Volterra  series  can  be  used  to  characterize  nonlinear  sys¬ 
tems  in  which  the  present  output  depends  on  the  present  and  past  values 
of  the  input  function.  Such  systems  are  referred  to  as  "systems  with 
memory". 

(d)  The  kernel  functions  themselves  are  not  necessarily  independent 
of  each  other,  but  they  are  independent  of  the  system  input.  Thus,  a 
parallel  system  realization  of  the  form  shown  in  Fig.  2-2  describes  the 
nonlinear  system  being  studied.  For  obtaining  an  nth  order  response,  it 
follows  that  the  first  n  kernel  functions,  or  their  transform  domain 
description,  must  be  determined. 

(e)  Some  important  advantages  of  the  Volterra  series  method  are 
that  it  places  the  input-output  relation  in  explicit  form  and  allows  us 
to  think  of  the  system  in  terms  of  functional  blocks,  as  shown  in  Fi<j. 
2-2. 


Fig.  2-2.  Volterra  series  equivalent  system  representation. 

This  often  allows  us  to  deduce  qualitative  properties  about  the  system 
which  may  not  be  clear  from  the  differential/integral  equations  used  to 
describe  the  system. 

We  now  proceed  with  an  important  aspect  of  the  Volterra  series 
method:  its  applicability.  Much  of  the  recent  work  on  this  topic  was 
done  by  Brilliant  C4],  who  introduced  a  new  topological  space  to  define 
the  concept  of  a  continuous  system  mathematically  and  showed  that  such 
systems  could  be  approximated  by  the  functional  representation. 

Once  again  consider  the  nonlinear  system  of  Fig.  2-1.  Let  x^(t) 
and  x^^t)  be  the  inputs  to  the  system  which  produce  y^(t)  and  y2<t)  as 
outputs,  respectively.  We  define  the  distance  between  the  input  func¬ 
tions  and  the  output  functions  as  follows: 
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(2-5) 


d1<x1,x2)  =  |  f  Cx^  (t)  -  x2(t)DcIt|  ,  t1  >  0 
t-t1 


and 


d2<y1'y2)  =  "  y2(t>*  (2~6> 

where  y^(t)  =  SCx^(t)3  and  y2(t)  =  SCx2(t)3.  Intuitively,  we  can  say 
that  S  is  a  continuous  system  if  d^  is  small  when  d^  is  small.  Brilli¬ 
ant  C43  gave  the  following  mathematical  definition  of  a  continuous  sys¬ 
tem:  For  a  time-invariant  system  S  with  bounded  inputs  x^  ( t >  and  x^Ct), 
S  is  continuous  if  for  any  e  >  0,  there  exists  a  6  >  0,  T  >  0  (T  suffi¬ 
ciently  large,  6  sufficiently  small)  such  that  if  d^  <  4  for  0  £  t^  <_  T, 
then  62  <  e.  It  has  been  shown  C4D  that  if  S  is  continuous,  according 
to  the  above  definition,  then  for  any  e  >  0  there  is  a  polynomial  system 
S£  such  that,  for  a  bounded  input  x(t),  |SCx(t)3  -  SeCx(t)3|  <  e. 

It  should  be  noted  that  the  distance  functions  defined  in  eqns. 
(2-5)  and  (2-6)  are  not  unique,  and  that  other  functions  could  have  been 
chosen.  The  definition  of  continuity  stated  above  is  just  one  defini¬ 
tion  and,  therefore,  the  aforementioned  condition  for  S  to  admit  a  func¬ 
tional  series  representation  is  only  a  sufficient  condition,  and  not  a 
necessary  one. 

The  conditions  on  the  system  S  that  admits  a  Volterra  functional 
series  description  can  be  qualitatively  summarized  as  follows: 

(a)  Output  is  a  single  valued  function  of  the  input. 

(b)  S  is  time-invariant.* 

*The  Volterra  functional  series  has  been  applied  to  time-varying 
systems  also.  See,  for  example,  C39D. 

18 


(c)  Small  perturbations  in  the  input  do  not  produce  abrupt  changes 


in  the  output. 

2-3 .  Derivation  of  the  Volterra  Series:  An  Example 

The  basic  problem  in  applying  the  Volterra  series  method  to  non¬ 
linear  systems  is  the  determination  of  the  kernel  functions.  In  this 
section  we  present  an  approach  for  deriving  these  kernel  functions.  A 
direct  and  more  algorithmic  approach,  which  also  shows  the  similarity 
between  the  Volterra  series  method  and  the  perturbation  method  C173,  is 
given  in  section  3. 

We  consider  a  system  with  a  single  nonlinear  element  and  the  fol¬ 
lowing  differential  equation  which  provides  the  input-output  relation¬ 
ship: 


L.Cy(t)3  +  eUCEa  yn(t)3  =  x(t> 
1  2  n=2  n 


(2-7) 


where  and  L ^  are  linear  operators,  defined  as  follows: 


,m  .m-1 

L.  CO  =  -2— -  CO  +  b.  3  +  ...  +  bflCO 

1  dtffl  ™  1  dtm  1  0 


(2-8) 


and 


j  t 

L2CO  =  +  62/£Odt  +  c3tO  (2-9) 

“CD 

where  e..  and  b..  are  constants.  Notice  that  if  e1  =  e2  *  0,  we  have  a 
memoryless  nonlinearity  such  as  a  resistor;  when  or  e2  are  non-zero, 
we  have  a  memory-type  nonlinearity  such  as  a  capacitor  or  an  inductor. 


Thus,  eqn.  (2-7)  describes  a  large  class  of  nonlinear  systems.  We  now 
proceed  to  determine  the  various  kernel  functions. 

If  e  =  0,  eqn.  (2-7)  is  the  linear  differential  equation: 

L1 Cy (t)D  =  x (t )  (2-10) 

Denoting  the  impulse  response  of  this  differential  equation  by  h(t), 
then  eqn.  (2-10)  is  known  to  have  a  zero  input  solution  plus  a  zero 
state  solution.  That  is, 

t 

y (t)  =  zero  input  solution  +  Jh(t-T)x(t)dT  (2-11) 

t0 

and,  second,  the  steady-state  solution 

t 

y  (t)  =  Jh(t-r>x(T)dT  (2-12) 

“CO 

which  arises  as  a  limiting  case  of  eqn.  (2-11)  when  all  transients  have 
decayed.  The  forcing  function  x(t)  is,  in  this  case,  assumed  to  act  at 
all  times  t,  -«  <  t  <  «. 

We  now  derive  the  Volterra  series  for  eqn.  (2-7).  We  assume  that 
for  t  <  0,  x(t)  =  y(t)  =  0.  Equation  (2-7)  can  then  be  converted  into  a 
nonlinear  integral  equation  by  applying  the  inverse  operator  to  L^, 
which,  under  the  assumed  initial  conditions,  is  just  the  following: 

—1  00 

CL.(p)3  x(t)  *  fh(t-T)x(x)dT  (2-13) 

0 

where  p  =  -gp  and  h(-r)  =  0  for  t  <  0.  The  resulting  nonlinear  integral 
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equation  is: 


JN  ® 

h(t-t)Li£a  y^Ct)]  =/ h(t-T)x(T)dx  (2-14) 

0  2  n=2  n  0 

We  now  assume  that  y(t)  is  given  by: 

y(t)  -  y^(t>  +  y£(t)  +  ...  +  y^(t)  +  (2-15) 

where  y^(t)  denotes  an  i-th  order  output.  Recalling  that  the  order  of 
an  output  is  defined  as  the  number  of  times  the  values  of  the  input  is 
multiplied  by  itself  in  the  convolution  integral,  we  immediately  get: 

y.(t)  =/h(t-x)x(T)dT  (2-16) 

1  0 

which  is  the  linear  system  response.  The  second  order  output  will  be 
due  to  the  squaring  of  the  first  order  output.  After  substituting  eqn. 
(2-15)  in  eqn.  (2-14),  we  get  the  following  integral  equation  involving 
the  second-order  output: 

y2(t>  +  e^h<t_t>L2i:a2y1(T)3dT  =  0  (2_17) 

Substituting  for  y^(t)  in  eqn.  (2-17)  from  eqn.  (2-16),  we  get 


yg^t)  =  “  e / h(t— 


(2-18) 


x(t^  JxCijJdi^dT^DdT 

cococo 

-  -  ea,  rf/hCt-TihCT-T^  )fi(T-T,)  • 

2  000  1  2 

x(t^ )x(T2)dt,|dT2dT  (2-19) 

where  h  is  obtained  after  the  Linear  operator  has  been  applied  to  the 
quantities  within  the  parenthesis  in  eqn.  (2-18).  When  =  0  in 

eqn,  (2-9),  i.e.,  for  the  memoryless  nonlinear  case,  fi  =  h. 

Similarly,  when  we  consider  the  third-order  terms  in  eqn.  (2-14) 
(which  will  be  due  to  the  second  order  and  first  order  outputs),  we  get 
the  following  integral: 

<®  _ 

y3(t)  +  e  Jh(t-x)L2Ca3y1  (t)  +  2  a2y-j<x)  y2(t)]dt  =  0  (2-20) 

"OP 

Substituting  for  y^(t)  and  y2(t)  in  eqn.  (2-20)  and  using  eqns.  (2-16) 
and  (2-19),  we  get 

y3<t)  3  “  eJh(t-T)Ca3J|0'fi(T“T^ )R(t"t2)R(t~t3) 


+  2a2  e/,0,/^T"’T1^^T"’a^^a”T2^ 


R(o_T3)do]x(T^)x(T2)x(T2)dT,jdT2dT2dT  (2-21 ) 

Higher  order  outputs  can  be  derived  in  this  bootstrapping  manner.  We 
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now  give  the  form  of  the  kernel  functions  for  the  nonlinear  system 
described  by  eqn.  12-7). 

When  n=1  in  eqn.  (2-4),  we  obtain  ^(t^)  from  eqn.  (2-16)  as: 

h^)  =  h(t)  *  i'1!——!  (2-22) 

-1 

where  i  denotes  the  inverse  laplace  transform.  It  is  noted  that  (s) 
is  the  characteristic  equation  of  the  linear  system. 

Setting  n=2  in  eqn.  (2-4)  and  comparing*  it  with  eqn.  (2-19),  we 
get 

<9 

h^l-r-j,^)  *  ~  ea2  J^h(r)R(T>j“t)fi(T2“T)dT  (2-23) 

Similarly,  setting  n=3  in  eqn.  (2-4)  and  comparing  with  eqn. 
(2-21),  we  get: 

QO 

^13^t1'12/13^  =  ”  erh(T)ea3fKTrt)fi(x2-T>fi(T3-r> 


?  • 

+  ^a2ej^Tl”T^^ff"’T^^T2~0^^T3"’a^c*a'^T  (2-24) 


As  mentioned  previously,  the  kernel  functions  are  independent  of 
the  input.  An  important  feature  of  the  kernel  functions  given  in  eqn. 
(2-23)  and  (2-24)  is  that  an  n-th  order  kernel  can  be  expressed  as  a 


*ln  eqn.  (2-1$T,  we  first  let  u  *  t-r,  u^ 
therefore,  =  t-u.,  t2  *  t-u2,  r*t.  * 

x~x2  *  u2“u*  *hen  «i  *  ti'  u2  s  x2' 

These  steps  are  worked  out  in  the  Appendix. 


=  t-t.,  u2  *  t-r2;  and, 
(t-u)-vt-u.)  s  u«*u  and 
u  *  t  in1  eqn.  (2-23). 


product  of  lower  order  kernals.  This  product  structure  property*  holds 
for  many  practical  lumped  system  problems.  For  example,  nonlinear  cir¬ 
cuits  with  polynomial  type  nonlinearities  have  such  kernel  functions. 
We  will  use  this  property  to  investigate  the  convergence  properties  of 
Volterra  series  for  this  important  class  of  nonlinear  system  problems. 

2-4.  Convergence  Properties  of  Volterra  Series 

The  previous  section  dealt  with  the  question  of  finding  the  kernel 
functions  from  a  nonlinear  differential  equation.  The  Volterra  series 
consists  of  an  infinite  sum  of  terms.  Clearly,  in  a  practical  applica¬ 
tion,  this  series  must  be  truncated  at  some  finite  term.  A  very  natural 
question  to  ask  is:  How  many  terms  in  the  series  are  required  to  give  a 
good  approximation  of  the  system  response?  Thus,  we  must  look  into  the 
convergence  properties  of  Volterra  series. 

The  question  of  convergence  of  Volterra  series  was  recently  exam¬ 
ined  by  Gilbert  C203,  and  Lesiah  C213.  Other  references  on  this  subject 
include  Barrett  C6,233,  Geyer  C34D,  and  Ku  and  Wolf  C223.  In  keeping 
with  the  application-oriented  spirit  of  this  dissertation,  we  first 
present  the  basic  convergence  theorem  from  C223  and  then  derive  some  new 
and  more  meaningful  parameters,  related  to  the  linearized  system,  that 
can  be  used  to  obtain  the  approximate  error  in  the  truncated  series 
solution. 


*Many  authors  £8]  refer  to  this  as  the  "separable"  property  of 
higher  order  kernels.  We  use  the  term  "product  structure"  instead, 
so  as  to  avoid  confusng  with  separable  kernel  in  linear  system 
theory,  which  denotes  a  different  concept. 


2k 
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a  a  |"l 

IjT***/  lhn<T1'*”'Tn)  *  I x Ct_Tp>  * dxp 

as  ao 

<  Xnf...f  |h  (t.,...,t  ) |dx. ...dx  (2-30) 

”  0  o  n  1  n  1  n 

Using  eqns.  (2-28)  and  (2-30)  with  (2-29)  gives 

|yCt>  I  <  I  Eyn(t)|  <  £snxn  (2-31) 

This  completes  the  proof. 

Equation  (2-31)  clearly  emphasizes  the  role  which  the  input  ampli¬ 
tude  plays  in  the  determination  of  the  number  of  output  terms  that  must 
be  retained  to  achieve  a  desired  accuracy.  The  importance  of  the  input 
amplitude  in  the  analysis  and  synthesis  of  nonlinear  can  never  be  over¬ 
emphasized.  The  existence  of  multiple  solutions  in  nonlinear  systems  is 
not  rare.  Whether  the  system  will  exhibit  the  desired  and  stable  solu¬ 
tion  very  often  depends  on  the  input  amplitude.  Thus  the  input  ampli¬ 
tude  plays  an  important  role  in  the  nonlinear  system  response.  Later  in 
this  section  we  shall  briefly  discuss  the  stability  and  boundedness  of 
solutions  via  the  Volterra  series  method. 

Returning  to  the  topic  of  convergence,  it  should  be  clear  that  eqn. 
(2-31)  is  not  of  much  use  in  its  present  form;  for  it  requires  the 
determination  of  the  n-th  order  kernel  before  the  bound  on  the  n-th  ord¬ 
er  response  can  be  found.  Even  when  the  kernel  is  determined,  the 
evaluation  of  Gn  may  not  be  a  simple  task.  In  an  engineering  problem 
one  prefers  to  have  an  a  priori  estimate  of  the  order  of  the  error  that 
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is  incurred  as  a  result  of  truncating  a  series  at  some  finite  term.  In 
the  following  paragraphs  of  this  sub-section,  we  explore  the  convergence 
of  lumped  systems  with  a  single  nonlinear  element  that  is  characterized 
by  the  differential  equation  (2-7).  We  then  investigate  the  convergence 
of  lumped  systems  with  multiple  nonlinear  elements. 

2-4.1.  Single  Nonlinear  Element  Case; 

In  this  section  we  study  a  nonlinear  differential  equation  with 
memoryless  nonlinearity  and  determine  the  bounds  on  the  solution  in 
terms  of  the  norm  of  the  linear  kernel  function,  and  concomitantly 
in  terms  of  the  poles  and  residues  of  the  linearized  system.  The  abili¬ 
ty  to  examine  the  convergence  of  the  solution  series  in  terms  of  the 
pole-residue  values  is  of  great  help  to  the  system  analyst,  because  such 
information  for  the  linearized  system  can  be  easily  obtained  from  many 
computer-aided  analysis  packages.  Furthermore,  the  results  derived  in 
this  section  give  a  means  of  estimating  the  bound  on  the  next  higher 
term  of  the  truncated  series  -  something  similar  to  the  case  of  numeri¬ 
cal  integration  formulae,  which  give  the  order  of  the  truncation  error 
in  terms  of  the  step  size. 

The  equation  to  be  examined  is  for  a  nonlinear  system  with  memory¬ 
less  nonlinearity  with  quadratic  and  cubic  terms.  We  examine  the  fol- 

* 

lowing  differential  equation: 

LCy (t)D  +  a2yZ(t)  +  a3y3(t>  =  x(t>  (2-32) 

where  L  is  a  linear  operator,  similar  to  in  eqn.  (2-8).  Then  using 
eqns.  (2-16),  (2-19),  and  (2-21),  we  can  immediately  write  the  approxi¬ 
mate  solution  for  y(t)  as: 
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y  (t)  =  y*h  (t~t)x  ( t)d  t  +  a2</jj[/h(t)h(T.j 


-t)h(T2_T) 


x(t^  )x(T2)dx^dT2clT  + 


2a^fJJJh(T>h(T1-T)h(a-T> 

0 


h  <  <  Tj~a)x  ( )x(x2)x(f3)dx,|dx2dx3dadx  + 

a2/jJfh(x)h(x,|“x)h(x2"T)h(xj-x)x(x^ JxCxg^xCx^)  • 


dT1dT2dT3dT  +  *** 

Now,  if  the  following  conditions  hold: 
|x(t)|  <  X  for  all  t 


(2-33) 


(2-34) 


and  / |h(x)|dx  =  G;  h(x)  =  0,  x  <  0  (2-35) 

0 

then  the  bound  on  eqn.  (2-33)  is  given  by: 

|y (t) |  =  Y  <  GX  +  |a2|G3X2  +  2|a2|265X3 

+  |a3|G4X3  +  •••  (2-36) 

In  fact,  if  higher  order  terms  were  considered,  the  bound  on  Y  will  be 
given  by: 


-m 
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(2-37) 


<  gx 

Y2  <  |a2|G3X2 

Y3  <  <2|a2|2G5  +  |a3|G4)X3 

\  <  CC4|a2|3  +  |a2|2)G?  +  C2|a2l |a3l  +  3|a3|)G6  +  |ajG53X4 


Y.  <  GC  E 
1  “  k=2 


ak|Yi,k] 


/ 


i  >  2 


(2-38) 


where 


Y 

i,  m 


i-m+1 

£  Y.Y  .  Y  =  Y 

j  —  y  j  t 


(2-39) 


Equations  (2-38)  and  (2-39)  give  a  recursive  relationship  for 
determining  the  bound  on  the  output  terms.  As  shown  by  the  first  few 
terms,  these  bounds  ultimately  depend  on  the  boundedness  of  the  norm 
of  the  kernel  function,  h^(i),  of  the  linearized  system  and  the  input, 
x(t). 

We  now  look  at  an  approach  to  estimate  G.  For  eqn.  (2-32),  the 
transform  domain  description  of  an  asymptotically  stable  first  order 
kernel  can  be  written  as*: 
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R. 


(2-40) 


♦Here  we  consider  the  case  of  linear  system  with  simple  poles.  The 
linear  system,  which  is  asymptotically  stable,  but  has  multiple  ord¬ 
er  poles  can  be  handled  similarly. 
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(2-44) 


and 


n/u . 

/  |g4<T)|dT 

0 


1  -  e 


o.t 

- 7 - ;  9^t)  =  R-e  sin.T 

iraj/aij  .1  J  J 


(2-48) 


Clearly,  6^  is  the  contribution  to  G  due  to  the  real  poles,  and  is 

the  contribution  due  to  the  complex  poles  in  the  linearized  system. 

Our  discussion  in  this  sub-section  has  dealt  with  the  convergence 
of  the  Volterra  series  for  nonlinear  lumped  systems  with  a  single 
memoryless  nonlinearity.  The  results  will  be  summarized  in  the  follow¬ 
ing  theorem,  whose  proof  is  included  in  the  appendix. 

Theorem  2-2:  Consider  the  equation 


L(^-)  •  y(t)  +  E2  a.y^t)  =  L.J  (—)  •  x(t> 


where  L(p)  =  pm  +  b  .pm  ^  +  •••  +  bn 

m- 1  U 

L1<p)  =  P1  +  CHPH  ♦••.•  +  c0  ,  A  <  m 

and  L(p)  is  assumed  to  have  roots  with  negative  real  part.  If,  for 
t  <  tg,  x(t)  =  y(t)  =  0,  and  for  t  21  tg,  |x(t)|  £  X,  then  y(t),  for 
t  21  tg,  remains  bounded  by  the  inequality  (2-37). 

It  should  be  noted  that  Theorem  2-2  shows  how  the  series  converges; 
it  does  not  show  if  the  series  converges. 

The  bound  on  the  output,  according  to  eqn.  (2-37),  ultimately 
depends  on  the  bound  on  the  1 ^  norm  of  the  linear  kernel  function.  This 
can  be  easily  obtained  by  the  method  discussed  above  for  finding  G.  If 
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the  series  converges,  these  resutts  can  be  used  to  determine  the  number 
of  terms  in  the  Volterra  series  that  must  be  retained  to  get  a  solution 
within  a  prescribed  accuracy;  one  does  not  have  to  first  determine  the' 
kernel  functions  (except  the  linear  kernel),  as  per  Theorem  2-1,  to  es¬ 
timate  the  number  of  terms  that  must  be  kept.  Knowing  6  and  X  and  the 
coefficients  a..,  the  bound  on  each  order  output  can  be  determined,  which 
tacitly  provides  us  the  information  about  how  many  kernels  need  be 
determined  to  meet  a  prescribed  error  criterion. 

Our  foregoing  discussion  has  dealt  with  nonlinear  systems  with 

memoryless  nonlinearity.  The  case  of  memory  type  nonlinearity  can  be 

treated  in  a  parallel  manner,  and  an  inequality  similar  to  eqn.  (2-37) 

derived.  The  only  difference  in  the  new  inequality  will  be  that  the 

norm  of  the  differentiated  or  integrated  first-order  kernel  function 

«0  , 

will  be  involved.  Thus,  for  example,  terms  of  form  f  |^Hi(t)|dT  or 


/  | h (x) | dx 
0 


dt  will  now  be  involved  in  determining  the  bounds  on  the 
various  order  responses.  These  terms  can  again  be  evaluated  easily. 
For  example,  consider  the  case  of  a  system  with  simple,  real  poles. 
Then, 


-p.t 

h (t)  =  Z  r.e  1  ,  p.  >  0 


i  1 


From  (2-49),  we  get 


S  =/  ( x)  I  d  t  -  |  h(x>  I 


(2-49) 


0 


We  now  consider  the  case  of  lumped  systems  with  multiple  nonlinear¬ 


ities. 

2-4 . 2  Multiple  Nonlinear  Element  Case 

The  case  of  lumped  systems  with  multiple  nonlinear  elements  is 
similar  to  the  case  of  one  nonlinear  element.  Instead  of  working  with  a 
single  differential  equation/  we  work  with  a  system  of  differential 
equations.  In  studying  the  convergence  of  the  Volterra  series  we  again 
use  the  "product  structure"  property  (PS)  of  the  higher  order  kernel 
functions  to  derive  the  bounds  on  the  various  order  outputs  in  terms  of 
the  norm  of  the  linear  kernels.  Before  proceeding  with  our  discus¬ 
sion  of  the  main  topic  of  this  sub-section  we  point  out  that  most  higher 
order  kernels  encountered  in  nonlinear  systems  with  polynomial  non- 
linearities  have  the  PS  property*.  In  some  cases,  a  "special"  viewpoint 
may  be  required  to  see  this  property. 

Consider,  for  example,  the  system  of  Fig.  2-3.  Van  Trees  C8D  uses 
this  as  an  example  of  a  system  with  non-PS  property  kernels.  His  claim 
follows  from  looking  at  the  second  order  input-output  relationship, 
which  is: 

z<tW  {  "H/VW  ’  Sb'VVV1  • 

x<t-T,j)x(t-T2>dT1dT2  (2-51) 


*The  PS  property  for  cascaded  systems  is  quite  apparent  from  the 
general  expressions  •  for  their  kernel  function.  See,  for  example, 
C53  to  verify  this. 
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Fig.  2-3.  A  System  with  "non-ps"  kernels. 


However,  if  we  view  the  output  of  each  nonlinear  element  separately,  and 
write  the  system  output  in  terms  of  these  component  outputs,  then  we 
find  that  the  kernels  of  the  system  of  Fig.  2-3  indeed  have  the  PS  pro¬ 
perty.  This  is  evident  from  eqn.  (2-51)  when  we  re-write  it  as: 

00  00 

=/  f  h*  ('Oh.  (x-)x(t“T1  )x(t-T,)dT.dT- 
0  0  1  a  1  1  a  2  1  2  12 

+  Hlb(x1)h1b(T2)x(t-Tl)x(t-T2)dT1dt2 

=  z1 (t)  +  z2(t)  (2-52) 

Then,  the  systems  relating  x(t)  to  z^(t)  and  x(t)  to  z2(t)  have  the  PS 
property.  Furthermore,  the  bound,  Z,  on  z<t>  is  given  by: 
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Z  =  |z(t)|  <  | Z 2 1  +  U2I 


(2-53) 


where  Z^  and  Z 2  are  the  bounds  on  the  component  outputs,  2^(t)  and 
z2(t),  respectively. 

We  now  study  the  convergence  of  lumped  nonlinear  system  by  working 
with  a  specific  example.  Consider  the  system  described  by  the  following 
set  of  equations: 

y\  *  y1  +  y1  +  y2  +  yz  +  a2yi  +  a3yi  =  x(t> 


*  2  3 

y*2  +  1 2  +  y2  +  y1  +  a4y2  +  a5y2  =  ^  (2-54) 

Just  as  for  the  case  of  single  equation,  we  apply  the  inverse  operator 
to  eqns.  (2-54)  to  obtain  the  integral  form  of  the  vector  equations: 


y1  (t)' 

f 

(t-x) 

h12(t-x>' 

a2y1^T^  + 

y2(t) 

h22(t-T) 

h22(t-x) 

a4y2(x>  +  a5y2(x) 

{ 


(t~ t) x (x) 

(t-x)x(i) 


dx 


(2-55) 


Now  we  assume  that  y^(t)  and  y-,(t)  are  made  up  of  various  order  outputs; 
that  is. 


y1<t)  *  y 
y2(t>  *  y 


(1) 

1 

(1) 

2 


where  y 


denotes  the 


(t)  +  y^Ct)  +  y^(t)  +  ... 
(t>  +  y22>(t)  +  y23)(t)  +  ... 

j-th  order  output  y^(t). 
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(2-56) 

Substituting  eqn. 


(2-56)  in  eqn.  (2-55),  and  collecting  terms  of  like  order  on  both  sides, 
the  first  two  terms  are  found  to  be: 


h.. (t~ t ) x (t )dx 
11 

(t“ t ) x  (t )dx 


we  can  readily  determine  the  bounds  on  the  first  and  second  order  out¬ 
puts,  provided  the  input  x(t)  is  bounded,  as  follows: 


(2-60) 


I y 2  ~  ^2  —  ^21^ 


r,  (2), ,  ,i 

L(2)l 

r  2  v2i 

|y1  (t)| 

Y1 

9nx 

i  (2) ,  , , 

= 

y(2) 

<  6 

2  „2 

|y2  «t)| 

1 

922X 

(2-61 ) 


The  above  example  of  a  system  with  multiple  nonlinear  elements  was 
deliberately  chosen  to  be  simple,  and  yet  we  find  that  the  determination 
of  the  bounds  on  the  output  is  quite  involved.  This  example  neverthe¬ 
less  illustrates  how  the  bounds  on  the  various  order  outputs  are  related 
to  the  matrix  involving  the  norm  of  the  linear  kernel  functions.  We 
summarize  the  results  by  the  following  theorem. 


Theorem  2-3:  Consider  the  system: 


L(^-)  £(t>  +  F t > U  =  x.(t> 


where 


L(p>  = 


*11<p)  a12<P>‘  *'•  a1n(p) 


*n1<p)  1n2<p>  •**,'nn<p> 


y(t>  =  £y<J(t>  y  g  (t  >  ...  yn<t)JT 
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'o'  >8 


=  <t> |  *  bound  on  the  t~th  order  output,  y  (t> 


) 


i  — t h  element  of  v_  is  1  if 

v  =  M  1  ...  0  ...  ID;  x^(t)  is  non-zero;  otherwise, 

it  is  zero. 

Proof:  The  proof  for  the  above  theorem  is  similar  to  that  for  Theorem 
2-2,  which  is  included  in  the  appendix. 

Theorem  2-3  gives  recursive  relationship  for  determining  the  bounds 
on  the  output  series  of  a  class  of  lumped  nonlinear  systems.  As  shown 
by  the  example,  prior  to  the  theorem,  these  bounds  ultimately  depend  on 
the  integrable  property  of  the  linear  kernel  functions,  h„(x).  An  es¬ 
timate  of  the  bound  of  these  kernel  functions  can  be  obtained  by  the 
method  outlined  in  section  2-4.1.  The  theorem  is  useful  in  determining 
the  highest  order  of  kernel  functions  that  are  needed  to  obtain  the 
solution  within  a  prescribed  accuracy.  Only  a  knowledge  of  the  G  matrix 
and  the  bound  on  the  input,  X,  are  required  for  this  task. 

It  should  be  pointed  out  that  the  basic  idea  in  the  above  theorem 
can  be  extended  to  nonlinear  systems  with  multiple  independent  and 
dependent  (coupled),  memoryless  or  memory-type,  nonlinear  elements. 
However,  the  notation  and  algebra  involved  becomes  so  complicated  that 
it  is  difficult  to  put  it  as  a  general  theorem.  One  can  nevertheless 
apply  the  idea  of  estimating  the  higher  order  responses  in  terms  of  the 
bounds  on  the  norm  of  the  linear  kernel  functions  and  the  inputs  for 
specific  system  problems  at  hand. 


2-4.3.  Stability  and  Boundedness  of  Solutions. 

Much  of  the  foregoing  discussion  in  this  section  has  dealt  with  the 
question  of  how  the  Volterra  series  solution  converges,  but  has  ignored 
the  question  of  whether,  if  at  all,  the  solution  for  a  system  converges. 
The  question  of  the  uniqueness  of  the  solution  has  also  not  been  ad¬ 
dressed  thus  far.  The  answers  to  these  later  questions  are  directly  re¬ 
lated  to  the  stability  and  boundedness  of  the  solution.  By  determining 
the  condition  on  X,  in  terms  of  the  ^  norm  of  the  linear  kernel  func¬ 
tion,  for  which  the  Volterra  series  converges,  we  can  get  much  insight 
into  the  stability  and  boundedness  of  a  solution.  Our  discussion  in 
this  sub-section  will  look  into  the  determination  of  the  constraint  on  X 
to  get  a  converging  series.  Theorem  2-2  will  be  most  useful  for  such  a 
task. 

Much  of  the  theory  on  the  uniqueness  and  boundedness  of  solutions 
of  ordinary  nonlinear  differential  equations  has  been  covered  in 
[15,21,353.  In  [153  the  nonlinear  differential  equation  is  converted  to 
the  Volterra  integral  equation  (which  is  an  equivalent  representation  in 
integral  form  of  eqn.  (2-7)),  and  using  the  Picard  iteration  approach, 
conditions  on  the  nonlinearity,  the  linear  kernel  function,  and  the 
boundedness  of  the  input,  are  derived  to  prove  convergence  and  unique¬ 
ness.  It  is  shown  that  if  the  forcing  function  is  bounded,  the  linear 
kernel  is  square  integrable,  and  the  nonlinearity  satisfies  the  Lip- 
schitz  condition,  then  starting  from  any  set  of  initial  conditions  at 
tg,  (2-7)  has  a  unique  solution  for  any  finite  range  of  time  t.  The 
proof  fails  if  either  tg  ♦  or  t  ♦  •.  Since  it  has  been  shown  C193 
that  the  sequence  of  iterates,  obtained  from  the  iteration  process. 
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correspond  to  the  partial  sums  o f  the  Volterra  series,  it  follows  that 
the  above  criterion  for  uniqueness  and  boundedness  can  also  be  applied 
when  the  Volterra  series  method  is  used.  However,  if  one  is  using  the 
Volterra  series  method,  it  is  only  natural  to  have  some  criterion  which 
establishes  a  direct  relationship  between  the  convergence  of  the  series 
and  the  stability  and  boundedness  of  the  solutions.  Most  authors 
120,213  have  established  conditions,  similar  to  those  stated  above,  to 
prove  convergence  and  uniqueness  of  the  Volterra  series  in  a  finite  time 
interval  C0,T];  the  case  of  CO,®],  or  C-®,®]  has  not  been  addressed. 
The  latter  problem  has  not  been  adequately  dealt  with  in  the  literature 
and  is  in  fact  a  difficult  problem.  The  only  known  work  in  this  area 
was  done  by  Barrett  C23].  We  use  results  from  C23]  to  a  specific  equa¬ 
tion  to  show  the  possibility  of  determining  the  stability  and  bounded¬ 
ness  via  the  Volterra  series  method. 

Consider  the  following  second-order  nonlinear  differential  equation 
with  a  cubic  nonlinearity: 

y  +  a^y  +  atf  +  ey^  =  x(t)  (2-61) 

Then  using  theorem  2-2,  we  can  immediately  write  the  bound  on  the  solu¬ 
tion  y(t)  as: 

|y (t)  |  =  Y  <  SX  +  MgV  +  ...  (2-62) 


where 


*♦3 


(2-63) 


G  =  / |h(T)|dt  ;  h(t)  =  t 
0 

X  =  max  | x(t)  | 

0<t<® 

Re-writing  (2-62)  as: 

Y  =  f (X )  =  GX  +  |e|G4X3  +  ...  (2-64) 

It  should  be  evident  that  if  the  series  for  f(X)  in  eqn.  (2-64)  con¬ 

verges  for  some  X  _<  X£,  then  the  Volterra  series  solution  for  eqn. 

(2-61)  also  converges  for  max  |x(t)|  <  X  ,  since  f(X)  dominates  the 

0«tO  c 

Volterra  series  term  by  term.  We  now  seek  a  nonlinear  algebraic  equa¬ 
tion  for  which  the  series  f(X)  is  a  solution.  The  roots  of  this  alge¬ 

braic  equation  along  with  the  constraints  on  X  for  which  the  series  is 
convergent  will  give  more  insight  into  the  solution  of  eqn.  (2-61).  All 
this  will  become  evident  as  we  proceed  along. 

It  can  be  shown  that  the  nonlinear  algebraic  equation*  for  which 
f (X)  is  a  solution  is: 

Y  -  |ej6Y3  =  GX  (2-65) 

Being  a  cubic,  eqn.  (2-65)  will  have  three  solutions  for  Y  in  terms  of 
X.  This  shown  in  Fig.  2-4,  where 


s  +  a,s  +  aQ 


♦After  the  inverse  operator  is  applied  to  eqn.  (2-61),  there  results 
a  striking  resemblance  with  eqn.  (2-65). 
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(2-66) 


Fig.  2-4.  Roots  of  the  Nonlinear  Algebraic  Equation 


Figure  2-4  provides  some  interesting  insight  into  the  behavior  of  the 
solutions  of  eqn.  (2-65),  and  concomitantly  into  the  solutions  of  eqn. 
(2-61).  First,  when  X  ♦  0,  then  Y  =  GX  -  which  implies  a  "linear  sys¬ 
tem"  operation.  When  X  <  X  ,  all  three  solutions  of  eqn.  (2-65)  are 

C 

given  by*: 


(2-68) 


The  range  of  convergence  for  the  solution  Y  in  eqn.  (2-68)  is  that 
-1  _ 

*If  sin  is  given  its  principal  value,  the  relation  Y  *  f(X) 
results. 
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| X |  <  X  ,  where  X  is  defined  in  terms  of  6  and  e.  Since  the  algebraic 

C  C 

equation  (2-64)  dominates  the  Volterra  series  solution  to  eqn.  (2-61) 
term  by  term,  it  follows  that  the  series  for  eqn.  (2-61)  will  also  con¬ 
verge  provided  | x (t ) |  <  X^,  for  all  t,  and  also  remain  bounded  by  the 
inequality  (2-62). 

The  above  example  illustrates  an  interesting  method  for  arriving  at 
the  constraint  on  the  input  function  in  order  to  insure  that  the  Volter¬ 
ra  series  solution  converges.  Instead  of  working  with  the  functional 
■series,  we  work  with  a  nonlinear  algebraic  equation.  The  behavior  of 
the  roots  of  the  algebraic  equation  in  a  given  region  then  provides  in¬ 
sight  into  the  convergence  of  Volterra  series. 

The  above  approach  will  also  show  the  cases  in  which  the  Volterra 
series  method  will  not  yield  meaningful  results.  Although  not  proven 
yet,  it  appears  that  if  the  linearized  system  is  not  asymptotically 
stable,  then  the  Volterra  series  solution  may  not  converge.  This  con¬ 
jecture  can  be  seen  when  we  consider  the  effect  of  a  simple  pole  which 
is  very  close  to  the  jw-axis.  This  pole  will  predominantly  effect  the 
^  norm  of  the  linear  kernel  function.  If  X£  has  G  in  the  denominator, 
as  in  the  above  illustration,  it  should  be  clear  that  as  G  +  •,  X^  ♦  0. 
Thus,  the  bound  on  the  input  for  which  the  series  converges  will  go  to 
zero  as  the  pole  moves  very  close  to  the  jw-axis. 
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CHAPTER  3 


MULTIPLE-NODE,  MULTIPLE-NONLINEARITY  CIRCUIT  ANALYSIS 

3-1 .  Introduction 

Nonlinear  systems  that  admit  a  Volterra  series  description  are  com¬ 
pletely  characterized  by  their  nonlinear  impulse  response  functions  or 
the  generalized  transfer  functions,  which  are  the  multi-dimensional 
transforms  of  the  nonlinear  impulse  response  functions.  Thus,  any 
analysis  of  nonlinear  systems  via  the  Volterra  series  method  will  entail 
the  determination  of  either  one  of  these  functions. 

Previous  works  C7,10,12, 13, 14,433  for  determining  the  generalized 
transfer  functions  have  relied  on  the  "harmonic  input"  method  or  the 
"nonlinear  current  source"  method.  The  harmonic  input  method  involves 
the  use  of  exponential  functions  as  inputs,  equating  terms  of  like  order 
on  both  sides  of  the  system  equation(s),  and  then  determining  the 
transfer  functions.  The  nonlinear  current  source  method  involves  the 
solving  of  the  nonlinear  differential  equation(s)  by  repeatedly  solving 
a  linear  differential  equation  with  nonlinear  excitation,  again  assuming 
harmonic  inputs.  Both  these  methods  give  a  recursive  relationship 
between  an  nth  order  transfer  function  and  up  to  the  (n-D-order 
transfer  functions.  Although  both  methods  are  mathematically  correct 
and  yield  the  correct  generalized  transfer  functions  for  a  system,  they 
both  have  a  common  "drawback":  there  is  no  direct  relationship  between 
the  approach  used  and  the  definition  of  the  generalized  transfer  func¬ 
tions. 
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Since  the  generalized  transfer  functions  are  transforms  of  the  im¬ 
pulse  response  functions,  it  only  stands  to  reason  that  the  determina¬ 
tion  of  these  transfer  functions  should  somehow  entail  the  use  of  an  im¬ 
pulse  function  as  an  input  to  the  system  -  just  as  in  the  case  of  linear 
systems.  Based  on  this  intuitive  approach  and  multi-dimensional 
transform  theory,  we  develop  a  new  method  for  the  determination  of  gen¬ 
eralized  transfer  functions  which  does  not  involve  the  use  of  exponen¬ 
tial  functions  as  inputs. 

As  one  would  expect,  the  results  obtained  from  the  approach 
developed  here  coincide  with  those  obtained  previously.  However,  it  is 
believed  that  this  new  approach  has  some  advantages:  first,  the  approach 
relies  on  fundamental  definitions  of  multi-variable  transforms,  and  is 
therefore  more  rigorous  than  the  "heuristic"  frequency  probing  methods; 
secondly,  the  morass  of  algebra  involved  in  the  harmonic  input  method  is 
avoided;  third,  the  new  approach  is  more  intuitive  and  provides  a  direct 
link  to  the  work  done  elsewhere  C53  on  nonlinear  systems. 

In  this  section,  we  develop  this  new  approach  gradually.  In  sec¬ 
tion  3-2  the  multi-dimensional  transform  theory  is  introduced,  along 
with  the  application  of  the  theory  to  specific  examples  which  will  be 
subsequently  used  in  deriving  the  generalized  transfer  functions.  In 
section  3-3  the  generalized  transfer  functions  for  an  r-th  order  non¬ 
linear  differential  equation  are  obtained.  Section  3-4  is  devoted  to 
the  determination  of  the  nonlinear  transfer  functions  of  a  general 
multiple-node,  multiple-nonlinearity  circuit  with  a  single  input.  The 
case  of  multiple  input  sources  is  treated  in  section  3-5.  A  specific 
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nonlinear  circuit  example  is  presented  in  section  3-6 


(3-7) 


f  ( t  ,  1 2, « •  •  ,t  ^)  —  ^  C F ( s ^ , s 2, • • • ,s^) ] 

Whether  we  use  Fourier  transform  or  Laplace  transform  in  eqns. 
(3-3)  and  (3-4)  depends  on  the  contours  of  integration  and  values  of 
s.j,s2,. .  .,sn.  The  importance  of  the  region  of  convergence  when  dealing 
with  unstable  and  non-causal  linear  systems  is  well  known.  Here  we  as¬ 
sume  that  the  systems  under  consideration  are  causal;  that  is,  the  Vol- 

terra  kernels  h  (t.,t-,. . ,,t  )  =  0,  for  t.,t,,...,t  <  0.  Also,  in  gen- 
n  12  n  1  2  n  — 

eral,  we  are  concerned  with  functions  (or  generalized  functions)  whose 
region  of  convergence  includes  the  imaginary  axis  in  each  variable,  so 
that  the  Fourier  transform  is  included  in  our  definitions. 

It  should  also  be  noted  that  most  of  the  properties  of  the  one¬ 
dimensional  transform  (linear  case)  carry  over  to  the  multi-dimensional 
case.  The  validity  of  this  statement  can  be  checked  elsewhere  C53. 

It  is  often  desirable  to  express  the  multi-variable  function, 
f (t^t^,. ../t^),  as  a  simple  function  of  time,  f(t),  and  vice  versa.  If 
all  t^'s  are  restricted  to  be  identical  so  that  t  =  t1  =  t2  =  •••  =  tn, 
then  f (t.j,t2,...,tn)  becomes  f(t).-  Thus,  in  the  two  variable  case,  f(t) 
can  be  obtained  from  f(t^,t2)  by  evaluating  f(t^,t2)  along  the  45°  line 
t1  =  t2.  Similarly,  if  we  plot  f(t^,t2,t^)  in  a  three-dimensional 
space,  then,  to  obtain  f(t),  we  are  only  interested  in  f(t^,t2,tj)  along 
the  line  t^  *  t2  *  ty  The  idea  of  converting  a  nonlinear  function  of 
one  variable  t  into  a  product  of  linear  multi-variable  functions  will  be 
used  repeatedly  in  the  sequel.  One  must,  however,  bear  in  mind  that  the 
ultimate  goal  is  to  obtain  the  solution  of  the  differential  equation  as 
a  function  of  time,  t,  and  that  the  introduction  of  t^,t2,  etc.  are 
merely  for  mathematical  manipulations. 
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We  now  apply  multi-dimensional  transforms  to  some  specific  cases 


which  will  be  subsequently  used  in  sections  (3-3)  and  (3-4). 

3-2.1  Volterra  Series:  The  Volterra  series  relates  the  system  input  x(t) 
to  the  system  y(t)  as  follows*: 


y(t)  =  £/•*•/  h  (t  ,...,t  )  n  X(t-T.)dT. 

n=1n-fold  n  1  n  i=1  11 


(3-8) 


where 


n 

yn(t)  =  /•**/  h  (t,,...,T  )  n  X(t-T.)dT. 

n  S-fold  n  1  n  i=1  11 


(3-9) 


Introducing  dummy  variables  t^t^^.-.^t^  in  eqn.  (3-9)  we  can  write 


yn(t)  as: 


yn(t)  =  yn(tvt2'-*’'tn>,t1=t.=...=t  =t 

i  c  n 


=  J"*/  n  X(t.-T.)dT 

n-fold  n  1  n  i=1  1  1  1 


(3-10) 


Taking  the  n-dimensional  transforms  of  eqn.  (3-10),  we  get: 


♦Unless  otherwise  stated,  all  limits  of  integration  are  between  0 
and  •  in  our  discussion  here. 


Yn(s1'“,/Sn>  =  «.Cyn(t1^...#tn>3 


/  n  -s .  t . 

=  h  (T  ,X  )  n  x(t.-T  )e  1  ’dT.dt.  (3-11) 

2n-fold  n  1  2  n  i=1  1  1  11 

Defining  tn-rn  =  a n,  t^-T^  =  =  o  ,  and  therefore: 

t  =  o+T/t  .  =  a  „,...,t.  -  CJ.+T  • 

n  n  n'  n-1  n-1  n-1'  '  1  11' 

d<»n  =  dtn,  =  dt^,.  ..,d =  dt^.  Substituting  these  quantities 

in  eqn.  (3-11)  and  performing  the  2n-fold  integrations  with  respect  to 

dt.  and  da^,  we  get 


YnCs1',**'sn)  =  Hn(s1'***/Sn)i^1X<si) 


(3-12) 


and  therefore  the  transform  domain  description  of  eqn.  (3-8)  becomes: 


w  n 

^^i,S)/‘",sn)  —  S  H  ( s. , ■ . . ,s  )  n  X(s . ) 
12  n  n=1  n  1  n  i=1  1 


(3-13) 


where  X(s)  is  the  transform  of  x(t).  If  the  input  x(t)  is  a  delta  func¬ 
tion,  then  eqns.  (3-12)  and  (3-13)  reduce,  respectively,  to: 


Yn<s1'*”'sn}  =  Hn(Sl,s2,...,sn) 


(3-14) 


and 


00 

Y(s(|,...,sn)  =  I)^Hn(s1,...,sn) 


(3-15) 


Equations  (3-12)  through  (3-14)  will  be  used  repeatedly  in  section 
(3-3) . 
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3-2.2  Nonlinear  Terms.  The  characteristics  of  nonlinear  elements  en¬ 


countered  in  many  nonlinear  dynamical  systems  can  be  represented  over 
any  finite  range  by  a  polynomial.  This  gives  rise  to  nonlinear  dif¬ 
ferential  equations  with  polynomial  type  nonlinear  terms.  When  such 
elements  are  used  in  a  system,  the  equilibrium  equations  contain  in¬ 
tegrals  and  derivatives  of  the  polynomials.  In  this  section  we  apply 
multi-dimensional  transforms  to  these  nonlinear  terms  to  obtain  general 
forms  of  the  transform-domain  description  of  these  terms.  The  details 
of  some  of  these  derivations  are  given  in  Appendix  B. 

y2(t)  Term: 


y  (t>  =  ycti>y<t2>lt1-t2-t 


(3-16) 


Y(srs2)  =  tCy2(t)3  *J/y<t1>y(t2>*  1  1  2  2dt.,dt2 


=  Y(s1)Y(s2> 


(3-17) 


y^(t)  Term: 


yvt)  =  y(t1)y(t?)y(t-)|t  . 

i  *  5  3 


(3-18) 
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d  n/iN  _ 

(t)  Term: 


d  n,^  _  d 


dt*  (t)  =  dty(t1'*‘*'tn)lt1*t-=...=t  ;  y(t1 . V  *  nyV  (3'23> 

it  n  psi 


YCs1#s2,...,sn>  =  aC£j^(V‘**'tn)dti-1 

S=  I  s 


=  (s.  +  s,  +  ...s  )Y(s„)...Y(s  ) 
td  n  1  n 


(3-24) 


(t)dt  Term; 


/yn(t)dt  =  Jyn<Vt'Vt"‘-'Vt)dt  (3-25) 

Letting  x..  -  t  =  t.,  and  taking  the  transform  of  eqn.  (3-25),  we  get 

Y(s.,s-,...,s  )  =  Y  (s4,...,s  >/(s,+s_+...+s  )  (3-26) 

i  c  n  n  i  nic  n 

The  general  forms  in  eqns.  (3-21),  (3-24)  and  (3-27)  will  be  used  in 
sections  (3-3)  and  (3-4).  The  salient  feature  in  each  of  these  equa¬ 
tions  is  how  an  nth  degree  polynomial  function  in  the  time-domain  is 
represented  by  the  nth-order  product  of  the  transform  of  the  function  in 
the  transform  domain.  It  is  this  product  structure  which,  analogous  to 
the  case  of  linear  system  analysis,  makes  the  analysis  of  nonlinear  sys- 
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terns  easier  via  the  transform-domain  approach. 


3-3 .  Nonlinear  Differential  Equation: 

In  this  section,  we  present  a  method,  based  on  applying  the  multi¬ 
dimensional  transforms  to  nonlinear  differential  equations,  to  determine 
the  response  of  a  nonlinear  system  with  a  functional  power  series  type 
of  nonlinearity.  The  nonlinear  differential  equation  considered  is  the 
following: 


L1  Cy<t)D  +  Ll^a/a)]  =  x(t> 
n=2 


(3-28) 


where  x(t)  and  y(t)  are  system  input  and  output,  respectively,  L1  is  a 
linear  differential  operator: 

R  Ar 

L.CO  =  E^-rC-3  (3-29) 
1  r=0dtr 

and  l_2  is  f,  or  a  constant,  or  a  sum  of  these  operators.  It  should 
be  noted  that  the  linear  operator,  Lg,  operates  on  a  polynomial  function 
of  y(t) . 

We  now  present  an  approach  whereby  the  nonlinear  differential  equa¬ 
tion  (3-28)  is  solved  by  a  bootstrapping  operation  by  first  dissolving 
it  into  a  set  of  linear  differential  equations  with  nonlinear  inputs. 
Multidimensional  transforms  are  then  applied  to  these  new  equations  to 
obtain  the  Volterra  series  solution. 

There  are  many  different  methods  of  rendering  a  nonlinear  differen¬ 
tial  equation  into  a  set  of  linear  differential  equation  with  nonlinear 
inputs.  We  use  the  approach  outlined  in  C14,173. 
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Assume  that  the  input  in  eqn.  (3-28)  is  of  the  form 

x (t )  =  ev(t)  (3-30) 

The  dummy  variable  e  helps  to  keep  track  of  the  order  of  the  terms:  a 
term  with  coefficient  en  signifies  an  nth  order  term.  This  can  be  seen 
easily  by  substituting  eqn.  (3-30)  in  eqn.  (3-9),  which  yields: 

/n 

h  (t,,...,t  )  riv(t-T.)dT.  (3-31) 

n-fold  n  1  n  i=1 

Let  us  assume  that  r(t)  is  the  response  to  the  input  v(t)  in  eqn  (3-28). 
Then,  according  to  the  Volterra  series  expansion,  as  per  eqn.  (3-8)  and 
(3-9),  the  n-th  order  response  is: 

n 

r  (t)  =  /.../  h(x_,...,T)  n  v(t-T.)dt.  (3-32) 

n  n=fold  n  1  n  i=1  11 

Comparing  (3-32)  and  (3-31),  we  obtain  the  following  relationships: 

yn(t)  =  enrn(t)  (3-33) 

and  therefore,  as  per  eqn.  (3-8), 

00  00 

y(t)  =  Ly  (t)  =  £e°r  (t)  (3-34) 

n=1  n  n=1  n 

We  now  have  two  differential  equations  which  relate  r(t)  and  v(t). 
First,  equation  (3-28)  can  be  re-written  as: 
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(3-35) 


L.  Cr(t)3  +  L5C£a  rn(t)3  =  v(t) 
i  c  “■>  n 


second,  after  substituting  eqn.  (3-34)  into  (3-28),  we  get: 

<»  N  » 

L.C  £  enr  (t)3  +  L.C  £  (a.  £  eV  (t))J3  =  ev(t)  (3-36) 

1  n=1  n  2  j=2  3n=1  n 

Thus  in  order  to  solve  eqn.  (3-28),  we  can  solve  eqn.  (3-36)  for  rn(t), 
n  =  1,2,...  and  substitute  in  eqn.  (3-34)  to  solve  for  y(t)  after  set¬ 
ting  e  =  1.  Setting  c  =  1  implies  that  x(t)  =  v(t),  and  therefore  y(t) 

=  r(t)  =£r  (t).  The  introduction  of  e  is  a  mathematical  artifice  which 
n  n 

helps  to  equate  coefficients  of  en  on  both  sides  of  eqn.  (3-36),  thereby 
yielding  linear  differential  equations  (involving  successively  higher 
order  outputs)  with  nonlinear  inputs.  This  is  simitar  to  the  perturba¬ 
tion  method,  [17,293. 

To  solve  for  r^(t),  the  linear  system  response,  we  equate  coeffi¬ 
cients  of  e*  on  both  sides  of  eqn.  (3-36),  thus  yielding  the  following 
equation: 

L1Cr1 (t)3  =  v (t)  (3-37) 

2  3  4  5 

Similarly  we  equate  coefficients  ofe,e,e,e,  and  so  on,  on  both 
sides  of  eqn.  (3-36)  to  obtain  the  following  equations: 

l^lr2U)l  +  L2Ca2r2(t)3  =  0  (3-38) 

^ Cr3(t)3  +  L2C2a2r<|(t)r2Ct)  +  a^t)]  *  0  (3-39) 
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L1Cr4(t)D  +  L2Ca2(2r1(t)r3(t)  +  rij(t))  +  3a3r*(t)r2(t) 

+  a4r^(t)D  =  0  (3-40) 

L1Cr5<t>3  +  L2C2a2r1  (t)r4(t)  +  a3(3r^  (t)r3(t)  + 

+  3ri(t)r2(t))  +  4a4r^(t)r2(t)  +  8^(03  =  0  (3-41) 


To  solve  for  the  generalized  transfer  functions  of  eqn.  (3-35),  we  take 
the  1 -dimensional  transform  of  eqn.  (3-37)  and  obtain: 

L1(s1)R1(s1) .  =  V(s1)  (3-42) 

If  v(t)  =  «(t),  then  V(s.j)  =  1,  and  therefore,  according  to  eqn.  (3-14), 
we  have 

=  »,<*,>  =  r^r  <3-43> 

To  solve  for  the  second-order  transfer  functions,  H2(s^,s2),  we  perform 
a  2-dimensional  transform  of  eqn.  (3-38)  to  obtain 

(s,j+s2)R2(s.|,s2)  +  ajLg^+Sg)*.)  (s^  )R^  (s2>  =  0  (3-44) 

? 

Notice  that  we  have  used  eqn.  (3-21)  to  transform  the  r^(t)  term.  Using 
(3-14)  and  (3-43)  in  eqn.  (3-44),  we  obtain 
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(3-45) 


a-L-Cs.+s.JH. (s* )H. (s,) 
22121112 

LjTs^j^Tp 


Simitar ly,  talcing  a  3-dimensional  transform  of  eqn.  (3-39),  we  obtain 


L1  (s1+s2+s3)R3(s1's2'S3>  +  L2(s1+s2)C2a2Ri{si)R2(s2/S3) 


+  a3R1(s1)R1(s2)R1(s3)D  =  0 


(3-46) 


Again,  using  eqns.  (3-43),  (3-45),  and  (3-14),  we  get 


R3  ^s^Sg^Sj)  H3  (s^ ,s2,s3)  ~  -  L2 (s^ ^Sg'^Sj)  C2a2H^  (s^  )H2 (s2,s3) 


+  ajH^j  (s^  )H^  (s2)H1  (s3)3/L,|  (s^Sg+Sj) 


(3-47) 


In  a  similar  manner,  we  can  derive  by  inspection: 


^4 ^S1 /S2'S3'S4^  —  L2 (  ^  )  Ca2 (2H,|  (s^  )H3 (s2,s3,s^) 


+  H2(s.j,s2)H2(s3,s^))  +  3a^H^  (s1>H1  (s2>H2 (s3,s^> 


+  a,  II  H1(s.)]/L1(  Es.) 

i=1  11  1  i=1  1 


(3-48) 


and 


*5  ^si's2's3's4's5^  "  ”  ^2^ ,  ^  ^  ^*1  ^4^2,*3,^4/*5^ 
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+  3a3(H1 (s1>H1 (s2)H3(s3,s4/s5) 


+  <s.j>H2(s2/s3>)H2<s^,s5>  ) 


+  ^a^H1 (s^ )H^ (s2)H1 (s3>H2<s^/s5>+a5 


5  5 

n  H- <s.>a/L, <  £  s • ) (3-49) 
1=1  1  1  1  i=1  1 


The  overbars  in  eqn.  (3-47)  to  (3-49)  represents  symmetrization 
operation,  which  is  an  averaging  operation.  In  general,  the  generalized 
transfer  functions  are  not  symmetrical  in  their  arguments;  that  is, 
H2(s^,S2)  may  not  be  equal  to  ^(s^s^).  The  symmetrization  operation 
on  an  unsymmetrical  system  is  performed  by  summing  each  of  the  nth  order 
transfer  function  over  all  permutations  of  its  arguments  and  dividing  by 
the  number  of  components  in  the  sum. 

The  use  of  symmetric  transfer  functions  is  not  merely  for  notation- 
al  convenience,  but  is  necessitated  by  the  method  we  use  for  introducing 
the  parameters  t<j,t2,...,  before  taking  the  transforms.  To  illustrate 
this,  we  note  that  v^Ovgtt)  can  be  written  as  v^  (t^ ^(tg,^), 
v.|  (t2)v2<t^,t3),  or  v^  (t3)v2(t,j,t2>.  The  first  term  has  transform: 
V.j  (s.j  )V2(s2,s3);  the  second  term  has:  vi  <s2)V2Cs,,Sj);  and  the  third  has 
transform:  V1  (s^IV^s^Sg).  When  V2(*,*)  is  not  symmetrical  in  its  ar¬ 
guments,  each  transformed  quantity  above  will  yield  a  different  value. 
Thus,  it  becomes  necessary  to  use  symmetric  transfer  functions  when  per¬ 
forming  numerical  computations  to  obtain  the  system  response.  It  can  be 
shown  that  the  response  is  unchanged  when  symmetrized  transfer  functions 
are  used.  For  example,  expressing 
v.j<t)v2<t)  *  j  Cv^  (t1)v2(t2,tJ)  +  v1(t2)v2(t1,tj)  +  v1(tj)v2(t.j,t2)3. 
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and  recalling  that  t^  =  t£  =  t3  =  t,  does  not  change  the  contribution 
due  to  v^lOvglt)  in  the  system  response.  In  the  remaining  part  of  this 
report  we  will  assume  the  generalized  transfer  functions  to  be  symmetric 
in  their  arguments. 

To  conclude  this  sub-section,  we  summarize  the  approach  for  obtain¬ 
ing  the  generalized  transfer  functions  of  a  nonlinear  system  and  also 
comment  on  the  important  ramification  of  the  method.  By  introducing  a 
dummy  variable  in  the  nonlinear  differential  equation  characterizing  the 
system,  a  set  of  differential  equations  of  the  following  form  was  ob¬ 
tained: 


LCrn(t)3  +  flr^lt))  =0,  n  =  2,3,...  (3-50) 

where  L  is  the  linear  system  operator  and  f(*>  is  a  nonlinear  function 

of  rn-1(t),  rn_2(t),...,r1 <t).  r^(t)  is  the  first-order  response,  which 

is  simply  the  response  of  the  linear  system.  The  relationship  in  eqn. 

(3-50)  is  clearly  a  recursive  one,  and  can  be  used  to  solve  for  r  (t) 

n 

in  terms  of  rn_^(t),  rn_2^'  etc.  This  is  d°ne  b/  taking  the  n- 

dimensional  transform  of  eqn.  (3-50)  to  solve  for  Rn(s^,...,sn),  which 

is  identically  the  nth-order  transfer  function  when  the  input  v(t)  is  an 

impulse  function.  The  transform  of  f(»)  is  done  by  inspection  with  the 

help  of  the  results  of  section  (3-2).  The  n-dimensional  transform  of 

LCr  (t)3  is  shown  to  be  L(s.+s,+...+s  >R„(s.,s,,. . .,s  ).  With  all  this 
n  i  c  nnl  c  n 

information,  eqn.  (3-50)  is  easily  solved  for  the  generalized  transfer 
functions. 

Much  of  the  discussion  above  has  been  concerned  with  the 
transform-domain  description  of  the  nonlinear  system.  It  should  be 
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pointed  out  that  the  form  of  eqns.  <3-3 7)  to  (3-41)  is  quite  suitable 


for  obtaining  the  time  domain  kernels  of  the  system.  By  defining  an 
operator  we  can  easily  solve  for  h.,(T.,).  Knowing  h.,(T.,),  we  can 

successively  solve  for  ^(t,,^),  etc.,  yielding*: 

h1  (t,)  =  Kit.,)  (3-51) 

h2(i;trt2)  =  a2/K(T)h1(t1-i)h1<t2-T)dT  (3-52) 

*13^T't1't2't3^  =  /K<T)Cajh^  (t,-t)h,  (t2-t)h  .  (t^-r) 

+  2a2h1(t1-T)h2<T1-T;t2/t3)DdT  (3-53) 

h^ ( i/t, ,t2#t3/t^)  =  -  ^K ( t)  Ca2 12h,  (t, -i)h3  ( t,, -t; t2,t3,t4) 

+  h2<T2"t;t1't2)h2<T3~T;t3't4>>  + 

3a3h, (t^-x)h^ (t2-T)h2<T4-T;t3^t4) 


4 

+  a,  n  h1(t.-T)]dx  (3-54) 

i=1  1  1 


with 


K(t)  i  l”1C1/L1  (s)3  (3-55) 

The  fifth-order  kernel  can  also  be  written  by  inspection  of  eqn.  (3-41).  . 


*Here  we  have  assumed  l2  to  be  a  "constant”  multiplication  operator, 
thus  giving  a  polynomial  type  memoryless  nonlinearity. 
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However,  for  the  sake  of  space,  we  do  not  write  it  here.  We  note  here 
that  there  is  a  one-to-one  correspondence  of  terms  between  the  time- 
domain  and  the  transform-domain  description  of  the  kernels  as  presented 
above.  Thus,  knowing  the  transform-domain  functions,  the  time-domain 
kernels  can  be  obtained  directly  using  the  above  equations.  The  pro¬ 
cedure  must  start  from  the  first  order  kernel  and  proceed  successively 
upwards  to  determine  the  second-order,  third-order, ...,n-th  order  ker¬ 
nels. 

3-4.  Multiple-Node,  Multiple-Nonlinearity  Circuit  Analysis 

Many  analysis  and  design  problems  in  circuits  and  systems  involve 
one  or  at  most  a  few  nonlinear  elements  in  an  otherwise  linear  time- 
invariant  circuit  or  system.  When  a  single  nonlinear  element  is 
present,  the  differential  equation  (3-28)  and  the  material  of  section 
(3-2)  will  be  adequate  for  analyzing  the  nonlinear  circuit.  For,  in 
such  a  case,  the  linear  circuit  can  be  characterized  by  a  convolution 
kernel  (via  the  Thevenin  or  Norton  Theorems)  to  give  the  overall  Volter- 
ra  integral  equation  C153,  which  can  also  be  cast  in  a  differential 
equation  form,  similar  to  eqn.  (3-28). 

However,  when  multiple  nonlinear  elements  are  imbedded  in  an  other¬ 
wise  linear  time-invariant  circuit,  the  analysis  entails  the  solution  of 
a  system  of  nonlinear  differential  equations.  The  approach  developed  in 
section  (3-2)  for  the  scalar  case  is  still  applicable,  but  must  be  ex¬ 
tended  to  solve  the  system  of  nonlinear  differential  equations. 

The  number  of  equations  to  be  solved  depends  on  the  number  and  the 
type  of  nonlinear  elements  considered.  When  only  independent  type  non¬ 
linear  elements  are  considered,  the  number  of  equations  is  less  than  or 


equal  to  the  number  of  nonlinear  elements  (assuming  that  the  output  is 
across  one  of  the  nonlinear  elements;  otherwise,  an  extra  equation  re¬ 
lating  the  nonlinear  element  voltages  (currents)  and  the  output  voltage 
(current)  is  needed  to  solve  for  the  output).  The  nonlinear  differen¬ 
tial  equations  in  such  a  case  is  again  derived  by  obtaining  the  Thevenin 
(Norton)  equivalent  circuit  (for  the  linear  part  of  the  nonlinear  cir¬ 
cuit)  at  each  of  the  ports  at  which  the  nonlinear  elements  are  present. 
When  dependent  type  nonlinear  elements  are  also  allowed,  the  analysis 
becomes  more  complicated;  for,  in  such  a  case,  the  controlling  vari¬ 
ables,  which  may  be  across  a  linear  element,  must  be  solved  for  and  sub¬ 
stituted  in  the  differential  equation  for  the  nonlinear  element. 

Previous  works  C7,10,143  for  determining  the  generalized  voltage 
ratio  transfer  functions  of  lumped  nonlinear  circuits  have  applied  the 
harmonic  input  method,  mentioned  previously  in  section  3-2,  to  the  nodal 
analysis.  C133  introduces  the  application  of  harmonic  input  method  to 
hybrid  analysis.  The  latter  approach  is  more  general  in  that  it  allows 
for  both  current-  and  voltage-controlled  nonlinear  elements. 

Our  discussion  in  this  section  for  solving  multiple-node, 
multiple-nonlinearity  circuits  will  be  centered  around  the  application 
of  multi-dimensional  transforms  to  a  cutset  type  analysis.  Thus,  we 
will  be  solving  for  the  generalized  voltage  ratio  transfer  functions. 
As  we  proceed  with  our  discussion,  it  will  become  apparent  that  a  cutset 
analysis  approach  is  the  most  natural  way  of  solving  for  the  generalized 
voltage-ratio  transfer  functions.  We  now  develop  the  procedure. 

The  first  step  in  the  analysis  is  to  represent  each  nonlinear  ele¬ 
ment  by  a  polynomial  expansion.  Thus,  in  the  distortion  analysis  of 
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transistor  amplifiers  C7,363,  the  exponential  type  controlled  sources  in 
the  Ebers-Moll  model  are  first  represented  by  a  Taylor  series  expansion 
of  the  function  about  the  quiescent  point,  thereby  yielding  a  polynomial 
in  terms  of  the  incremental  variables.  The  types  of  nonlinear  elements, 
and  their  series  representation,  that  are  commonly  encountered  are: 

1.  No  memory,  independent  nonlinearity  (Nonlinear  Resistor) 

CO 

i  =  F(v)  *  Da.v3  (3-56) 

j=1  3 

2.  No  memory,  dependent  nonlinearity 

CO  00 

i  =  G(u,v)  =  £  Ea.Vv*  ,  ann  -  0  (3-57) 

j=0k=0  Jlc  DU 


3.  Capacitive,  independent  nonlinearity 


1  =3fa(''>  4s-i«i 

at  dt j=1  j 


A.  Inductive,  independent  nonlinearity 

i  *  J$(v)dt  =  J  £a.v3dt 
--j=1  3 

where 

i  =  incremental  current  through  the  element 
v  s  incremental  controlling  voltage 


(3-58) 


(3-59) 
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u  =  incremental  controlling  voltage 

The  general  procedure  employed  to  solve  for  the  nonlinear  transfer 
functions  of  a  single-input/  single-output  nonlinear  circuit  using  the 
cutset  analysis  approach  is  illustrated  in  Fig.  3-1  by  considering  each 
of  the  four  nonlinear  element  types  mentioned  above. 

Consider  the  nonlinear  circuit  N,  shown  in  Fig.  3-1 (a),  containing 
a  nonlinear  resistor,  a  nonlinear  dependent  source,  a  nonlinear  capaci¬ 
tor,  and  a  nonlinear  inductor,  where  each  nonlinear  element  is  voltage 
controlled.  The  procedure  begins  by  identifying  all  the  nonlinear  ele¬ 
ments,  as  shown  in  Fig.  3-1 <b).  We  note  that  the  four  nonlinear  ele¬ 
ments  depend  on  six  voltages.  The  next  step  is  to  lump  the  linear  parts 
of  the  nonlinear  elements  with  the  existing  linear  network  to  form  the 
augmented  linear  network.  The  square,  cubic,  quartic,  etc.  terms  of  the 
nonlinearity  are  treated  as  nonlinear  current  sources,  indicated  by  i£, 
meaning  the  nth  order  current  source  at  port  k.  Since  the  dependent 
source,  glv$,v6),  depends  on  voltages  v$  and  v^,  we  also  extract  these 
as  ports.  Thus,  altogether  we  end  up  with  an  8-port  linear  network,  as 
shown  in  Fig.  3-1 (c). 

The  output  variables  to  be  found  are  the  voltages  at  these  eight 
ports.  The  augmented  linear  network  is  denoted  by  N'  in  Fig.  3-1 (c). 
To  solve  for  the  voltage  vector  v_ =  Cv^  Vg  Vj  ...  Vg3,  we  immediately 
recognize  that  the  branches  across  these  voltage  variables  must  be 
selected  as  part  of  the  tree  C29],  Clearly,  some  of  the  other  branches 
in  the  augmented  linear  network  may  also  appear  as  part  of  the  tree. 
These  will  then  appear  as  voltage  variables  in  the  cutset  equations  for 
the  augmented  linear  network.  Since  there  is  no  need  for  these  addi- 
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tional  variables,  we  can  reduce  the  dimensionality  of  our  equations  by  a 
systematic  elimination  of  these  unwanted  variables.  In  the  case  under 
consideration,  we  should  be  left  with  only  the  vector  v.  =  Cv-j  v^  ...  Vg3 
as  the  unknown  vector.  Each  of  these  8  ports  will  have  a  set  of 
transfer  functions  of  order  1  to  n  associated  with  it.  Our  task  here  is 
to  solve  for  these  transfer  functions. 

At  this  point,  we  make  the  following  general  notational  defini¬ 
tions: 


H^(Si,S2^«»*/S|()  — 


Hk(Sl,...,sk> 

2 

Hk  ,  *  *  *  /S({  ^ 


(3-60) 


where 


Hj|  =  kth  order  nonlinear  transfer  function  from  the  input  to  the  jth 
port;  m  =  8  in  our  example  here. 


v(t>  *  Cv.(t)  v,(t)  ...  v  (t)j 
—  «  c  m 


(3-61 > 


where  v.  =  voltage  at  the  ith  port 


The  cutset  equations  for  the  m-port  nonlinear  network  can  be  writ¬ 


ten  as 


(3-62) 


Y(p)v^  +  F  (v)  +  Gtu^)  +  pQ.lv)  +  ~  ^(v)  = 

tv  /z  (p)DC1  0  0  ...  o:T 

9  9 

where 

p  =  differential  operator, 

Y(p)  =  Reduced  admittance  matrix  for  the  p-port  augmented  linear 
network 

F_(v_)  =  vector  composed  of  all  nonlinear  currents  through  the  zero 

memory  independent  nonlinearity 

G(u,v)  =  vector  composed  of  alt  nonlinear  currents  through  the  zero 
memory  dependent  nonlinearities 

GKv)  =  vector  composed  of  all  nonlinear  currents  through  the  non¬ 

linear  capacitive  nonlinearities 

4<v>  =  vector  composed  of  all  nonlinear  currents  through  the  non¬ 

linear  inductive  elements. 

z  (p)  =  source  impedance 

•J 

Since  the  linear  parts  of  the  functions  F(»),  G(»),  Q(»),  and  ♦<•)  in 

eqn.  (3-56)  through  (3-59)  have  been,  lumped  together  with  the  linear 

part  of  the  network,  the  general  form  of  these  functions  will  be  as  fol¬ 
lows: 
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•  •  • 


(3-63) 


Z/v)  =  2^  ( v )  +  Z_3  ( v )  +  Z^(v)  + 


where 


ZgCv)  is  a  quadratic  function  of  v_ 

Z^(v)  is  a  cubic  function  of  v_ 

Z^(v)  is  a  quartic  function  of 

•  •  • 

MO  being  MO,  £lO,  0.(0,  or  $(•).  Thus,  eqn.  (3-62)  can  be  re¬ 
written  as: 


Y (p)v  =  — 
~  —  z 


V  (t) 

a 

o 

o 


’  lk(t) 


,  k  >  2 


(3-64) 


where  j_^(t)  denotes  vectors  of  2nd  and  higher  order  current  sources  due 

1 

to  Liv2s  and  ^(v).  T^*e  mathematical  artifice  used  in 

section  (3-2)  could  have  been  applied  here  also  to  obtain  the  form  of 
all  the  nonlinear  current  source  terms,  j_^(t).  For  the  sake  of  brevity, 
we  will  not  use  that  approach  here,  but  simply  use  the  results  of  sec¬ 
tion  (3-2)  to  identify  the  different  order  current  sources  due  to  dif¬ 
ferent  nonlinearities.  These  are  summarized  in  Table  3-1,  where  v1(t) 
denotes  the  ith  order  response  voltage  v(t),  which  control  the  nonlinear 
element  characteristics. 
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Table  3-1.  Nonlinear  Current  Sources  in  multiple-node,  multiple- 
nonlinearity  circuit  analysis. 

Nonlinear  Resistor,  : 

k  =  2:  a2Cv^32 

k  =  3:  2a2Cv1v2]+a3Cv1]3 

k  =  4:  a2C2v1v3+(v2)23+3a3Cv1:2v2+a4i:v1D4 

Nonlinear  Dependent  Nonlinearity  : 

.  _  -  r  1 ,2  r  1  _2  11 
k  =  2:  a20Cu  3  +  aQ2Cv  ^  +  a^u  v 

k  =  3:  a3gCu^  33  +  agjCv1!3  +  a^Cu^v*  +  a^u'cv'*]2  +  2a2Qu1u2  + 

,  1  2X  r  1  2  2  1. 

2aQ2v  v  ♦  a^Cu  v  +  u  v  3 

-  /  r  In4  x  r  In4  x  V  1n3  x  r  Vr  V 

k  =  4:  a^gLu  3  +  ag^Cv  3  +  a-j3u  *-v  3  +  a22^u  ^  ^v  3  + 

,,  13,  P  2n2.  .  ,31,  13.  2  2,  ,  ..,13.,. 

a^  (2u  u  +  Lu3)  +  a^(u  v  +  u  v  +  u  v  )  +  ag2(2v  v  + 

r  2-,2,  x  T  r  I-,2  2  x  ,  r  1-I2  2  x  ,r  I,2  2  x 

Cv  3  )  +  3a3gCu  3  v  +  3bq3Lv  3  v  +  a^  (Cu  3  v  + 

,  1  2  1,  .  ,  2r  1-2  .  ,  1  1  2. 

2u  u  v  )  +  a^2<u  Cv  3  +  2u  v  v  ) 

Nonlinear  Capacitive  Nonlinearity  £fl  ( v> : 

,  -  r  1n2 
k  =  2:  a2pCv  3 

k  *  3:  2a2pCv1v23+a3pCv133 

k  *  4:  a2p(2v1v3+Cv232)+3a3pCv^32v2+a4pCv^34 


I 


73 


Table  3-1  (contd.) 


Nonlinear  Inductive  Nonlinearity,  C1_/£^#(v) 


2  12 

k  *  2:  — ^CvV 

P 


.  _  2a2,  1  Z.Jlr  S3 

k  =  3:  — -iv  v  J+-— Cv  1 
P  P 


k  =  4:  ^i<2v1v3+Cv2]2)+^£v132v2-t^Cv1D4 
P  p  p 


We  observe  that  the  nonlinear  current  source  terms  in-  Table  3-1  are 


similar  to  the  nonlinear  terms  whose  transforms  were  derived  in  section 
3-2,  except  for  the  nonlinear  dependent  source  terms,  which  are  func¬ 
tions  of  two  controlling  voltages  u  and  v.  The  form  of  the  transforms 
of  the  nonlinear  dependent  source  will,  however,  be  similar  to  the  other 
nonlinearity  types.  These  can  again  be  written  by  inspection.  For  ex¬ 
ample, 

a2otu1(t)]2  =  a2ou1 <t ^ lu1 (t *20  u^si)u(s2)  (3-65) 
a^u1  (t)v^  (t)  =  a^u^  (t.j )v^  (t2>  -M'  a^-j  Uls^VlSg)  (3-66) 
a20u1  (t)v2(t)  =  S2QU1  (t^  )v2(t2/tj)  «-*  a2Q  U(s1)V(s2,s3)  (3-67) 

and  so  on. 

We  also  note  that  a  k-th  order  current  source  term  in  Table  3-1 
depends  on  responses  of  order  less  than  k,  which  implies  that,  in  order 
to  calculate  a  transfer  function  of  order  k,  we  need  to  determine  the 
transfer  functions  up  to  order  (k-1). 

The  first  order  transfer  function  can  be  solved  for  easily.  It  is 
simply  the  linear  circuit  response.  Therefore, 

Y(p)v(t)  =  i)(t)  <3-68) 

For  a  single  input  system,  j.(t)  =  1/z  Cv  (t)  0  0  ...  03T,  where 
vg(t)  is  the  source  voltage.  Taking  the  transform  of  eqn.  (3-68),  and 
assuming  that  the  input  source  to  be  an  impulse  function,  we  get: 


(3-69) 


V^s.)  =  Ms.)  =  1/z„  CY  Cs^  )D_1  Cl  0  0  ...  0]T 

—  i  —ii  g  ~  1 

where  Kjls.j)  was  defined  in  eqn.  (3-60). 

The  equation  for  obtaining  the  second-order  response,  as  per  eqn. 
(3-64),  is  the  following: 

Y(p)v(2)(t)  =  -  j^( t)  (3-70) 

Since  the  input  to  the  nonlinear  circuit  is  assumed  to  be  an  impulse 

function,  the  transform  of  eqn.  (3-70),  after  using  eqn.  (3-14),  is: 

Y(s^+S2)H2^s-j/s2^  =  (3-71) 

The  elements  of  vector  I^(s^,S2)  can  be  obtained  by  performing  a  two- 
dimensional  transform  on  the  terms  associated  with  k  =  2  in  Table  3-1. 
This  operation,  as  indicated  earlier,  can  be  carried  out  by  inspection. 
Thus,  we  have 

^(s^Sg)  =  ~  CYls^+Sg)]  ^(s.j^)  (3-72) 

Likewise  we  can  solve  for  Jfjls^Sg,^).  In  general,  we  solve  for  the 
nth  order  transfer  function  using  eqn.  (3-73): 


-1, 


(s^ ,S2, •  •  •  ,s^)  =  CY_(  22  s^)3  ^ (s^ , • ■ . ,s^) 


(3-73) 


We  observe  a  striking  similarity  between  eqn.  (3-73)  and  the  equa¬ 
tions  for  nodal  or  cutset  analysis  encountered  in  linear  circuit 
analysis.  A  little  thought  would  show  that  the  process  of  solving  eqn. 
(3-73)  is  identical  to  solving  the  linear  circuit  in  Fig.  3-2.  We  have 
nonlinear  current  sources  as  inputs  to  the  augmented  linear  circuit.  A 
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Augmented  Linear 
Network,  N' 


j-i  ^ 


7) 


+  h: 


Transfer  Functions 


k-th  order  vector  of  transfer  functions  is  obtained  by  exciting  the 
linear  circuit  by  the  kth  order  current  sources.  Just  as  in  the  case  of 
linear  systems,  superposition  can  be  applied  here  when  a  particular  ord¬ 
er  response  is  determined  from  the  lower  order  responses.  That  is,  a 
k-th  order  response  can  be  obtained  by  applying  the  k-th  order  current 
sources  one-by-one  at  each  of  the  ports  and  then  summing  up  the 
responses.  It  is  important  to  note,  however,  that  the  complete 
responses  of  order  up  to  (k-1)  must  be  determined  before  we  can  obtain 
the  kth  order  response  by  superposition.  It  is  also  noted  that  the 
illustration  of  Fig.  3-2  is  for  pedagogic  purpose  and  that  the  nonlinear 
current  sources  are  not  physically  present  in  the  circuit  under  con¬ 
sideration. 

3-5.  Multiple  Input  Circuit  Analysis 

Much  of  the  foregoing  discussion  has  been  concerned  with  the 
analysis  of  nonlinear  circuits  with  single  inputs.  However,  many  appli¬ 
cations  of  practical  significance  in  nonlinear  circuit  analysis  have 
multiple  inputs.  For  example,  in-  a  receiver  system,  the  mixer  circuit 
has  two  inputs:  1)  the  message  signal,  and  2)  the  local  oscillator  sig¬ 
nal.  The  transmitter  again  has  nonlinear  circuits  with  multiple  inputs. 
The  Volterra  series  method  is  especially  well  suited  for  the  analysis  of 
such  circuits.  In  this  section  we  discuss  how  the  various  order 
transfer  functions  change  as  a  result  of  multiple  inputs. 

From  the  discussion  in  section  3-4,  it  should  be  apparent  that  the 
analysis  of  nonlinear  circuits  using  the  Volterra  series  method  involves 
the  repeated  analysis  of  a  linearized  circuit.  The  fundamental  rela¬ 
tionship  had  the  following  form  (see  eqn.  3-64): 
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(3-74) 


Y(p)v^t)  =  2  ~]p^  j_-|  (t )  -  j_2<t)  -  Vj(t)  +  ... 

where  2^(t)  is  the  k-th  current  source  vector.  The  k-th  order  current 
source  (k  2),  as  per  our  previous  discussion,  depends  on  up  to  the 
(k-1)  order  voltage  ratio  transfer  functions,  is  injected  at  each  of  the 
ports  at  which  the  nonlinear  elements  are  present,  and  is  due  entirely 
to  the  nonlinear  characteristics  of  the  nonlinearity.  Furthermore,  it 
is  proportional  to  the  k  values  of  the  circuit  input  multiplied  togeth¬ 
er.  Thus,  the  number  of  elements  in  the  vector  ii^(t),  k  >_  2,  remain  un¬ 
changed  when  multiple  inputs  are  present;  only  the  v.(t)  vector  is 
changed. 

Consider,  for  example,  the  two-input  circuit  of  Fig.  3-3(a).  Then, 
to  solve  for  the  first-order  transfer  function,  we  write  the  vector 
transform  equations  as: 

y(s1)V(s1)  =  ^(s^  (3-75) 

where 

I1(s1>  =  CYg1(s1)Vg1(s1)  Yg2(s1)Vg2<si)  0  0]T  (3-76) 

and  Y  and  are  as  defined  previously.  The  transfer  function  vector  can 
be  written  as: 

H^s.,)  *  <H10Cs1>  +  H01(Sl>  (3-77) 

where 
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fig.  3-3.  (contd.)  Multiple  Input  Nonlinear  Circuit  Analysis 


— 10  <S1 5  = 


WC1 )  ,  .  ..(2).  . 

V  (s^)  V  (s1 ) 


gi 


gi 


v<p)(Sl) 

vgi 


(3-78) 


V  ,=0 

g2 


and 


fioi  (V  = 


v‘%)  v(2)(Sl) 


g2 


g2 


v(p)(Sl) 


gi 


(3-79) 


vgi=0 


where  is  the  voltage  at  port  i. 

The  second-  and  higher-order  transfer  function  vectors  are  solved 
for  by  removing  the  given  input  sources  and  aoplying  the  fictitious  non¬ 
linear  current  sources  across  the  ports  at  which  the  nonlinear  elements 
are  present.  The  vector  transform  equation  for  solving  for  the  second- 
order  transfer  function  is  still  given  by: 

H^s-^Sg)  =  -CY(s1  +  Sg)]*1  Cl2(srs2)D  (3-80) 

where 

l2(svs2)  =  Cl(1)(Sl,s2>  I(2>(srs2>  ...  I(p)(svs2)]  (3-81) 

(l) 

Depending  on  the  nonlinearity  type,  the  general  form  \>f  I  (s^,s2>,  the 
second-order  current  source  across  port  t,  will  be: 

I2t)(s1's2)  =  a2Hi<t)(si>Hia><s2>  (3-82) 

(t) 

where  (*)  is  known  from  eqn.  (3-77).  The  determination  of  the 
higher-order  transfer  functions  is  done  similarly. 

In  summary,  we  note  that  the  presence  of  multiple  input  sources  in 
a  nonlinear  circuit  does  not  drastically  alter  the  procedure  for  deter- 
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mining  the  Volterra  transfer  functions.  Only  the  structure  of  the 
first-order  current  source  vector  is  changed  as  a  result  of  multiple 
sources.  This  change  is  reflected  in  the  values  of  the  elements  making 
up  the  second-  and  higher-order  current  source  vectors,  whose  structure 
remains  unchanged. 

3-6.  An^  Example 

In  this  sub-section,  we  present  an  example  to  illustrate  the  ideas 
presented  in  the  previous  sub-sections.  Specifically,  we  derive  the 
nonlinear  transfer  functions  for  a  multiple  nonlinearity  circuit. 

Consider  the  circuit  of  Fig.  3-4(a),  which  contains  a  resistive  and 
a  capacitive  nonlinearity.  The  nonlinear  element  descriptions  are  also 
given  in  the  Figure.  The  augmented  linear  network  is  shown  in  Fig. 
3-4(b).  The  equilibrium  equations  for  the  nonlinear  circuit  can  be 
written  by  inspection  as: 


G+pC  -pC  0 

v^  (t) 

-pC  pC+pC^+1/pL  -1/pL 

v2(t> 

_°  -1/pL  1/pL+91_ 

v3Ct> 

Gv  <t)‘ 
s 

0 

0 


C2PV2  +  C3PV2 

2  3 

92v3  +  93v3 


(3-83) 


where  v.(t)  is  the  i-th  node  voltage,  p  =  and  G  =  1/R.  Each  of  the 
i  at' 

node  voltages  can  be  expanded  into  a  Volterra  series,  thus  giving 


(3-84) 


v.(t)  =  £  v{k)(t) 

k=1  1 

(k) 

where  (t)  represents  a  k-th  order  term  in  the  solution  for  the  i-th 
node  voltage.  Substituting  eqn.  (3-84)  in  eqn.  (3-83)  for  i  =  1,2,  and 
3,  and  equating  terms  of  like  order  on  both  sides  of  eqn.  (3-83),  we  get 
the  following  sets  of  equations  for  the  various  order  responses. 


First  Order: 


G+pC  -pC  0 

-pC  pC+pC^+1 /pL  -1/pL 

0  -1/pL  17PL+9l 


(IK  , 
(t) 

(IK  , 
v^  (t) 

(IK  N 
v3  (t) 


Gv  (t) 
s 

0 

0 


(3-85) 


Second  Order: 


G+pC  -pC  0 

-pC  pC+p^+1/pL  -1/pL 

0  -1/pL 


r  (2) ,  vi 

"  “ 

v1  (t) 

0 

«k  . 

v2  (t) 

1 

II 

c^Cv^’ct):2 

(2K  , 
v3  (t) 

(3-86) 


Third  Order: 


G+pC  -pC  0 

-pC  pC+pC^+1 /pL  -1/pL 


-1/pL  1/pL+9lJ 


v{3)(t) 

v«)(t)‘ 

v3a>Ct) 
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0 


2C2pCv 

2g2Cv 


(1)  (2) 
2  v2 


(1)  (2) 
3  v3 


x  „  r  (1  ),3 
3  +  CjpCv^  3 

,  r  (1  )n3 
3  +  g3Cv3  3 


The  higher  order  response  equations  can  be  written  similarly, 
solve  for  the  transfer  functions  from  the  above  equations. 

Taking  the  one-dimensional  transform  of  eqn.  (3-85),  we  get 


r  m  V*1 

G+s1c  "S1C  o 

V,  (.,> 

fGV  (sj“l 
s  1 

— s^ C  s,j  (C1+C)+1  /s^L  — 1  /s^L 

V21><s1> 

= 

0 

0  -1/s.jL  1  /s^L+g^ 

0 

Assuming  v  (t)  =  6(t),  and  therefore  V  (s.)  =  1,  we  get: 

S  Si 

H(1)(Sl)  =  CZ(s,j  )3C6  G  03T 

where 


(1)  V1<1>(s1>  V21>ls1)  V31>(s1>  T 

s  s  s 


and 


Z(Sl)  = 


G+s^ 


’S1C 


-s.C  s.  (C+C,)+1/s.L  -1/s.L 
1111  1 

0  1/s1L+g1 


-1 


To  solve  for  the  second-order  transfer  function 


,(2) 


H,  (s^/Sg),  we  must  first  recognize  the  following  with  (t^  *  tg 


(3-87) 


We  now 


(3-88) 


(3-89) 


(3-90) 


(3-91 ) 


vector, 
*  t): 
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00 

LCv{2>(t)D  =  L[//hf)(TrT2)Vs<t1-T1)Vs(t2-T2)dT1dT2] 

=  Hf2)(s1/s2>Vs(s1>Vs(s2)  (3-92) 

Since  vs^s^)  *  Vg(s2)  =  the  right-hand  side  of  eqn.  (3-92)  is  identi- 
(2) 

cally  H..  (s^,s2).  Similarly,  we  write  (with  =  t): 

Lkv”  >(t)32J  =  LC//vif1)(T1)vf1>(T2)«(t1-T1)«(t2-T2)dT1dT23 

=  vi<1)(*1)vj2>(*2)  (3-93) 

We  note  here  that  V?^(s4)  s  wf1)(s1)  since  v  (t)  =  fi(t). 

Till  s 

Taking  the  two-dimensional  transform  of  eqn.  (3-86),  and  substitut¬ 
ing  eqns.  (3-92)  and  (3-93)  in  it,  we  get: 

H<2>Cs1,#2)  =  -  Z(s1  +  s2>tO  (s1+s2)C2H21)(s1)H21)(s2) 

92H31)(s1)H31)<s2>]T  <3-94) 

where  H,  Z  are  defined  in  eqns.  (3-90)  and  (3-91). 

To  solve  for  the  third-order  transfer  function,  we  proceed  in  a 

similar  manner.  However,  as  mentioned  in  section  3-4,  symmetrization 
operation  is  used  in  solving  for  the  third-order  transfer  function  vec¬ 
tor.  The  symmetrization  operator  is  associated  with  the 
Cv^*(t)vj2*(t)3  term.  To  see  this,  we  write 
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rwC1),«v  (2) .  _  1  ,  (1),,.  .  (2K  .  .  .  (IK  .  <2K  .  . 

Lvi  (t)vi  (t)J  =  j  Lvi  (t^v.  (t2/t3>  +  v.  <t2)vi 


(3-95) 


trt2=t3=t 


2  vf1 1  <t1  >v^)  (t2^t3>  (3-96) 

We  note  that  when  we  let  t^  =  t2  =  tj  =  t  in  eqn.  (3-95),  the  two  sides 
are  identically  equal.  As  mentioned  previously,  the  variables  t..  's  are 
used  for  mathematical  manipulations;  the  final  analysis  involves  only 
one  variable  t.  Thus,  with  this  in  mind,  we  note  that  the  response,  as 
a  function  of  t,  will  remain  unchanged  when  symmetrization  is  performed. 

Proceeding  in  a  manner  similar  to  that  for  solving  for  the  second- 
order  transfer  function  vector,  we  obtain  the  following  third-order 
transfer  function  vector: 

(s^,s2,s3)  =  Z(s,|+s2+s3)  • 


(Sl+s2+s3) )  (Sl >H22)  (s2/s3)+C3H21  }  (s,  } <s2)H^  5 (s3) 


ag^1 }  (s^  )H32)  (s2,s3)+g3H31 )  Cs1 )«” }  (s^H*1  J  <Sj) 


(3-97) 


We  note  that  all  the  quantities  on  the  right  side  of  eqn.  (3-97)  are 
known,  and  thus  we  can  solve  for  the  third-order  transfer  functions. 
Our  derivations  here  have  been  fairly  detailed;  it  should,  however,  be 
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pointed  out  that  the  form  of  these  transfer  functions  can  be  obtained  in 
an  algorithmic  manner  by  inspection,  thus  rendering  Volterra  series 
method  a  viable  approach  for  computer  aided  analysis  of  nonlinear  cir¬ 
cuits. 
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CHAPTER  4 


STEADY-STATE  AND  TRANSIENT  ANALYSIS 


4-1 .  Int  roduction 

In  the  analysis  of  nonlinear  systems,  two  main  classes  of  solutions 
are  generally  sought:  1)  steady  state,  and  2)  transient  response.  In 
this  section  we  show  how  these  solutions  can  be  obtained  via  the  Volter- 
ra  series  method. 

In  this  section  we  present  how  the  steady  state  and  the  zero-state 
transient  responses  can  be  obtained  from  the  generalized  transfer  func¬ 
tions,  which  can  be  obtained  in  an  algorithmic  manner  -  as  was  shown  in 
section  3.  For  the  sinusoidal  steady-state  response  we  deal  with  the 
case  of  multi-tone  input;  for  the  transient  response,  we  treat  the  case 
of  input  signals  which  are  Laplace  transformable  and  factorable  -  which 
clearly  includes  a  large  class  of  signals  used  in  everyday  application. 
No  numerical  integration  is  required  in  obtaining  both  these  solutions. 

Section  4-2  deals  with  the  subject  of  obtaining  the  sinusoidal 
steady-state  solution  using  the  generalized  transfer  functions.  Section 
4-3  presents  the  approach  for  obtaining  the  zero-state  transient 
response. 

4-2.  Sinusoidal  Steady-State  Analysis 

In  linear  system  theory,  the  sinusoidal  steady-state  response  is 
intimately  tied  to  the  transfer  function  of  the  system.  A  similar 
result  is  found  for  higher  order  responses  using  the  Volterra  series 
method:  an  n-th  order  response  at  a  particular  frequency  is  directly  re¬ 
lated  to  the  n-th  order  transfer  function.  In  this  section  we  develop 
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this  relationship. 

If  the  harmonic  input  method  [7,10-143  had  been  used  in  deriving 
the  generalized  transfer  functions  in  section  3,  the  relationship 
between  the  n-th  order  steady  state  response  and  the  n-th  order  transfer 
function  would  have  been  self-evident.  But,  since  multi-dimensional 
transform  theory  was  used  to  derive  the  generalized  transfer  functions, 
this  relationship  must  be  developed.  We  treat  the  specific  case  of  n=2 
in  section  4-2.1  and  then  derive  the  general  relationship  in  section 
4-2.2. 

4-2.1 .  Second-order  Sinusoidal  response; 

The  second-order  output,  according  to  the  Volterra  series,  is  given 
by; 


y,(t)  =  //  h,(t-x  ,t-t  >x(t )x(t  )dT  dt. 
2  0  0  2  12  1  2  1  ‘ 


(4-1) 


Consider  the  input  signal  comprising  two  unit  sinusoidal  signals  at  fre¬ 
quencies  uig  and  The  input  x(t)  is  therefore; 


x(t) 


exp(j»T>+exp(-jttT>" 

exp  ( j  «bT )  +exp  (— j  &>b  t )" 

2 

h  2 

Substituting  eqn.  (4-2)  in  (4-1),  we  have; 


y?(t)  =  /  /  h-(t-t  ,t-x  )  • 
2  0  0  2  1  2 
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!®XP<j“aT1 
— 


)+exp(-j<u  Ti)  exp(j 

-7 - — 


jMbT1)+exp(-jo.bT1)1 

— ? - j 


exp(j“,aT2>+exP(“)t‘,aT2^  expljUbT^+expl-jWb-r^)' 
- 2 - - -  + - 2 - 


•  dT^d-Tg  (4-3) 

Considering  one  cross  term  only, 

•»  “  . 

/  J  h2(t-T1^t-T2)  exp(ja)aT1+jubT2>dT1dT2  (4-4) 

and  Letting  a ^  =  t-t^  and  a 2  =  t~r2  and  carrying  out  the  integration 
yields, 

H2(j<i)a,jwb)expCj  <(i»a+o)b)tD  (4-5) 

Considering  the  other  cross  term  similarly  yields 

-j- H2(j»b,j(fla)expCj  (ua+«b^^  (4-6) 

However,  if  H2(s,j,s2)  is  symmetrical  in  its  arguments,  as  they  are  as¬ 
sumed  to  be  in  this  dissertation,  then  the  terms  in  eqns.  (4-5)  and 
(4-6)  are  equal.  The  complex  conjugate  terms  appear  similarly.  Hence, 

the  output  at  frequency  «a+ulj  is: 

y  ( t )  |  ,  =  |  H—  ( j  a)  ,j«.  )|cosC(u  +0).  )t  +  6  .3  (4-7) 

<■>  +«.  2  a  b  a  b  a+b 

a  p 

The  2ufi  or  2ub  term  and  their  complex  conjugates  appear  only  once  in 
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eqn.  <4-3);  hence,  their 
*1 

j|H2(ju»b,ju)b)|,  respectively 
the  results  would  be  similar. 


magnitude  will  be  ,j«  )  I  and 

.  If  only  one  frequency  input  was  present. 
The  second-order  output  would  then  be: 


| HpC jw  ,ju  >| 

y2(t)  =  |H2<j(1)a,-jua)|  + — = - 2 - 2—  cos(2«at+82a>  (4-8) 

Thus,  if  we  know  H2(s^,s2>,  then  the  quantities  in  eqn.  (4-8)  can  be 
easily  evaluated.  This  is  analogous  to  the  case  of  linear  systems, 
where  the  complex  variable  s  is  replaced  by  ju  to  compute  the  response 
at  a). 

If  more  than  two-tones  were  present  at  the  input,  the  second  order 
response  would  be  evaluated  by  taking  all  combinations  of  two  frequen¬ 
cies  at  a  time. 

The  response  of  the  third  and  higher  orders  is  similarly  treated. 
We  now  present  the  general  case. 

4-2.2.  General  Sinusoidal  Steady-State  Analysis. 

In  this  sub-section,  we  develop  the  relationship  which  can  be  ap¬ 
plied  directly  to  compute  the  sinusoidal  steady-state  response  of  a  non¬ 
linear  system  from  its  nonlinear  transfer  functions,  which  can  be  ob¬ 
tained  by  the  method  presented  in  section  3.  The  discussion  here  relies 
heavily  on  C103. 

Consider  a  nonlinear  system  excited  by  the  sum  of  K  distinct  tones; 
i.e.,  defining  N  =  2K,  we  have. 


x  (t )  = 


_  1 


7  £  A.  exp(jw.t) 
c  i*1  1  1 


(4-9) 
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where  will  include  both  positive  and  negative  frequencies,  and  for 
a  negative  frequency  will  be  the  complex  conjugate  of  A.,  for  the  posi¬ 
tive  frequency  in  order  to  have  x(t)  real.  Then,  the  nth  order  output, 
yn(t),  is  given  by: 

y  (t)  =  /•••/  h  (X  ,T  )  n  x(t-T  )dT 
n  n-fold  n  1  n  i=1  11 


-/"•/ VV  —  'V  7  £  AkexpCj<uk(t-T.)]dT. 


(4-10) 


Carrying  out  the  product  operation  in  eqn.  (4-10),  we  get  a  function 
yn(t)  containing  Nn  terms,  given  by: 


N 

yn(t)  *  £ 

kri 


N  . 

E  -tA,  •••a.  h  ( jw, 


kn=1  2n  k1 


k  n 
n 


•  expCj  (a),  +...+U.  )t]  (4-11) 
1  Kn 

Notice  that  in  arriving  at  eqn.  (4-11),  we  have  performed  the  x..  in¬ 
tegration  in  eqn.  (4-10),  thus  giving  rise  to  the  n-th  order  transfer 
function  in  eqn.  (4-11).  As  the  indices  are  varied  over  the  range  1 
to  N,  many  of  the  terms  will  be  at  the  same  frequency.  The  number  of 
terms  at  various  particular  frequencies  will  vary  according  to  what  fre¬ 
quency  combinations  are  taken.  For  example,  in  the  case  of  n= 2  in  sec¬ 
tion  4-2.1,  there  were  two  cross  frequency  terms,  while  there  was  only 
one  second  harmonic  (at  2u>.)  term.  Similarly,  for  n*3,  there  are  six 
terms  in  eqn.  (4-11)  at  frequency  wa+u>b+wc,  three  terms  at  2ua+u>b,  one 
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term  at  3a>  ,  etc.  The  nonlinear  transfer  functions,  which  make  up  the 
coefficients  of  these  frequency  terms,  differ  only  in  their  arguments. 
However,  since  the  transfer  functions  are  assumed  to  be  symmetric,  the 
coefficient  of  the  output  at  frequency  w  +u>. +o»  (in  the  case  of  n=3) 

3  D  C 

can  be  multiplied  by  6.  This  obviates  the  need  for  taking  all  combina¬ 
tions  to  compute  the  output  at  <^+<1)^+0^.  Likewise  we  handle  the  case  of 
other  frequency  combinations.  With  this  insight,  we  can  peek  at  the 
problem  from  a  different  perspective. 

Let  m^n^,. .  .,mN  be  non-negative  integers.  Then,  the  number  of 
terms  at  frequency  =  miui+n,2a,2+***+lllV')N  is  equal  to  the  number  of 
ways  of  forming  m^u^+...+mNa)N.  In  the  n-th  order  output  spectrum  to  a 
multi-tone  input,  each  term  is  evaluated  by  taking  a  distinct  combina¬ 
tion  of  n  input  tones  at  a  time.  To  compute  the  n-th  order  output  when 
the  input  frequencies  are  “i'<1'2'**‘'“n/'  we  must  therefore  restrict  m..  in 
the  following  manner  to  compute  Wj.: 


"!  +  ra2  +  ’**  +  fflN  =  n 


(4-12) 


Now  the  problem  reduces  to  the  following:  find  the  number  of  ways  in 
which  n  objects  can  be  divided  into  N  groups  of  which  the  first  contains 
m^  objects,  the  second  m^  objects,  etc.  The  solution  to  this  problem  is 
given  by  the  multi-nomial  coefficient  C403: 


m.  !ra 


n! 


12* 


(4-13) 


By  deriving  eqn.  (4-13),  we  have  obviated  the  repetition  of  terms  that 
is  inherent  in  eqn.  (4-11).  An  equivalent  way  of  representing  eqn. 
(4-11)  through  the  use  of  eqn.  (4-13)  then  becomes: 


•  expCj(m1u1  +  ...  +  mNu>N)t3  (4-14) 

Since  yn(t)  is  real,  eqn.  (4-14)  also  contains  the  complex  conjugate 
terms.  Thus,  the  coefficient  of  the  sinusoidal  term  at  frequency 
m^iD,j+...+mNta)N  in  the  n-th  order  output  is  given  by: 


m.  m_ 

A1*2 
*1  *2 


rv~i  Hn(j<u1,...,j«2,...,. 


mutinies  mutinies 


m^times 


(4-15) 


In  computing  the  entire  n-th  order  response  in  eqn.  (4-14),  we  take,  all 
distinguishable  combinations  of  m,.  satisfying  eq.  (4-12).  According  to 
C40D  there  are 


n+N-1 


_  (n+N-1)! 
"  n! (N-1) ! 


(4-16) 


such  combinations. 

Equation  (4-14)  is  the  fundamental  relationship  between  the  n-th 

order  output  and  the  n-th  order  transfer  function.  At  first  glance,  the 

evaluation  of  this  equation  appears  to  be  a  formidable  task.  But,  after 

some  thought,  one  finds  that  this  is  not  such  a  difficult  task  after 

all.  We,  however,  defer  the  discussion*  of  this  till  section  5. 

*A  computer  implementation  for  n=2.3  appears  in  C413.  One  can  ex¬ 
tend  it  for  any  n,  with  the  ultimate  limitation  being  the  storage. 
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We  now  illustrate  the  use  of  eqn.  (4-15).  We  assume  that  the  non¬ 
linear  transfer  functions  are  known.  The  case  for  n=2  can  be  easily 
verified  from  the  discussion  in  section  4-2.1.  For  a  two-tone  input  at 
<0^  and  and  n=3,  we  have  the  following  cases: 

(a)  The  output  at  and  have  the  following  amplitudes,  respectively: 


3!|A2I‘ 


|  M3Cju»i  ^-j<i»2^  j<*>2>  I 


(4-17) 


3! A  |A  |2 

y3(t),a.2  -  amiTi'T  iH3<)“1^-j“1^“2>l 

(b)  The  output  at  2<»><|+u2  has  the  following  magnitude: 

31  A^a 


(c)  The  output  at  3«1  has  the  following  magnitude: 


>'3<t>l3.. 


3!  <An)' 
(4)3! 


(4-18) 


(4-19) 


(4-20) 


The  other  combinations  can  be  carried  out  similarly.  For  the  above 

cases  we  make  the  following  observations:  both  eqns.  (4-17)  and  (4-18) 

are  similar  to  obtaining  the  output  at  &>a+wb+wc,  and  therefore  we  see  a 

3!  (=6)  multiplication  factor*,  which  accounts  for  the  six  combinations 

at  <d„+w,  +»  that  were  mentioned  earlier;  eqn.  (4-19)  is  similar  to  ob- 
a  b  c 

*The  constant  factor  4  in  the  denominator  appears  consistently  in 
all  the  output  terms,  and  is  therefore  not  regarded  as  a  variable 
multiplication  factor  here.  This  factor  appears  due  to  the  way  x(t) 
was  expressed  in  eqn.  (4-9). 
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taining  the  output  at  2w  +qik ,  and  therefore  has  a  multiplication  factor 
of  3,  which  again  is  in  accordance  with  our  earlier  discussion;  eqn. 
(4-20)  is  like  evaluating  the  output  at  3u-,  and  hence  has  a  multiplica- 

a 

tion  factor  of  1. 

In  section  (3-5),  we  dealt  with  the  analysis  of  multiple  input  non¬ 
linear  circuits.  In  obtaining  the  sinusoidal  steady-state  response  of 
such  circuits  the  material  of  this  section  is  still  applicable.  Howev¬ 
er,  care  must  be  taken  in  keeping  track  of  the  various  input  frequen¬ 
cies,  and  their  associated  transfer  functions,  when  such  an  analysis  is 
warranted. 

4-3 .  Transient  Analysis  using  Volterra  Series 

In  linear  system  analysis,  the  most  fundamental  relationship 
between  the  input  and  the  output  is  the  convolution  integral.  Thus,  to 
compute  the  output  response,  one  merely  convolves  the  input  function 
with  the  impulse  response  function.  However,  when  the  input  function  is 
Laplace  transformable,  the  convolution  operation  is  rarely  performed  and 
one  resorts  to  the  simpler  transform  methods  for  obtaining  the  system 
response.  The  transform  method  involves  the  use  of  the  system  transfer 
function. 

A  scheme,  analogous  to  the  transform  domain  method  in  linear  system 
analysis,  can  be  devised  for  obtaining  the  response  of  a  nonlinear  sys¬ 
tem  via  the  Volterra  series  method.  Again,  the  inputs  allowable  under 
this  scheme  must  be  Laplace  transformable.  This  restriction  is  not  very 
severe  since  most  of  the  inputs  considered  in  everyday  applications  are 
Laplace  transformable.  The  scheme  presented  here  uses  the  generalized 
transfer  functions  and  the  transform-domain  description  of  the  input  to 
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compute  the  output.  The  response  thus  obtained  is  the  zero-state  time- 
response. 

The  basic  problem  considered  here  is:  Given  the  generalized, 
transfer  functions  of  the  nonlinear  system  and  the  transform  of  the  -in¬ 
put,  determine  the  output  response  without  performing  a  convolution  in 
the  transform  domain.  A  procedure,  due  to  George  C53,  helps  to  solve 
this  problem  by  associating  time  variables  in  the  transform  domain. 
That  is,  given  Y^s.j^)  as  the  transform  of  the  output  function 
y2(ti,t2>,  then  Y^Cs),  the  transform  of  y2<t),  will  be  found  directly 
from  Y2(s1's2>*  technique  is  called  "association  of  variables"  and 

is  applicable  to  the  class  of  lumped  systems.  Lumped  systems  have  the 
property  that  all  the  generalized  transfer  functions  are  factorable 
(just  like  the  linear  situation);  that  is,  if 

f(trt2>  =  f  1  <t1)f1  (t2>  (4-21) 

then 

F(srs2)  =  F(s1)F(s2)  (4-22) 

The  examples  in  section  3  -  in  particular  the  example  in  section  (3-6)  - 
illustrate  the  concept  of  factorable  transforms. 

In  the  following  sub-section  we  present  the  "association  of  vari¬ 
able"  technique.  In  section  4-3.2  we  apply  this  technique  to  the  time- 
domain  analysis  of  multiple  node,  multiple-nonlinearity  circuits. 


4-3.1.  Association  of  Transform  Variables 


The  association  of  variables  technique  151  is  derived  here  for  the 
second-order  response.  For  higher-order  transforms,  one  simply  associ¬ 
ates  variables  successively  using  the  results  for  the  second-order 
case. 

The  two-dimensional  transform  pair  is  related  as  follows: 


Y2^s1's2 


>  *  /  /  y2(t1,t2>exp(-s1t1-s2t2)dt1dt2 


C4-23) 


and 


y2(t1't2)  =  (2 7]")2//  Y2<s1/s2>exp(s1t1+s2t2>ds1ds2  (4-24) 

We  would  like  to  find 

y2(t1)  =  X2<t1,t1)  (4-25) 

without  actually  performing  the  inverse  transformation  indicated  by  eqn. 
(4-24) . 

The  transform  of  y2(t^)  is  Y2(s1),  where 

OB  St 

W  =  <dy)  /  Vs)e  1<js  (4"26> 

But  from  eqns.  (4-24)  and  (4-25)  we  have 

y2(V  =  y2(t1'V 
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=  (^T>2  //dslds2  WS2>e 


(4-27) 


(S1+S2H1 


Letting  =  s,  we  have 

W  a  sr  /  exp<st1>ds 


•  SJ  l  Vs’s2'S2>ds2 


Equating  (4-28)  and  (4-26) ,  we  have 


Y2<s)  '  2irj  /  Y2<s"s2's2>ds2 


(4-28/ 


(4-29) 


Equation  (4-29)  still  involves  a  convolution  operation,  and  hence  is  not 
too  useful  for  our  purpose.  In  the  case  of  lumped  systems,  where  the 
transfer  functions  are  factorable,  this  convolution  can  be  performed  by 
inspection.  For  example,  a  typical  second  order  transform  is  of  the 
following  form: 


Y  (s  s  )  =  *  •  1  _  •  t 

t2's1's2*  s^+Sg+x  s^+y  s2+z 

Then,  substituting  eqn.  (4-30)  in  eqn.  (4-29),  we  have: 

y2(s)  ~  ST  /  s-s2+s2+x  *  s-s2+y  *  s2+z  ds2 


(4-30) 


10) 


CO 

_  _x_  .  1  (  V 

S+X  F7  s-s2+y 


s,+z  ds2 


The  integral  in  eqn.  (4-31)  reduces  to: 

1  Y  Z ,  _  YZ 

“FT  -s^+7s+y>  s2+z  S2  ~  s+y+z 


Thus,  we  find  that: 


r/q;,  bec0""s 


and 


_Y _ 

s^y 


becomes 

s2+z 


Y2 

s+y+z 


(4-31) 


(4-32) 


A  list  of  commonly  encountered  transform  terms  in  lumped  nonlinear  sys¬ 
tems  along  with  their  transforms  after  the  association  of  variables  is 
given  in  Table  4-1.  A  more  detailed  list  can  be  found  elsewhere  C243. 


Table  4-1.  List  of-  Associated  Transforms 


Multi-dimensional  Transforms  Associated  Transforms 


F(s,+s,+...+s  )  F(s) 

12  n 

k  k 

(s^+a) (s2+b)  s+S+B 

k  k 


^s^+Sg+a)  (s,j+b>  (s2+c) 


(s+a)(s+b+c) 


(s^+Sg+Sj+a)  <s,j+b)  (s2+c)  (Sj+di 


Ts+aT7s+B+c+3T 
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Table  4-1.  <contd.> 

_ _  k  k  .  j_  r  i  .  ii 

(s^+s^+Sj+a)  (s^+Sj+b)  ^s^+ci  ^s^+c)  2a-b  s+a  |_s+b+c  s+3cj 

4^3._2  Transient  Analysis  of  Nonlinear  Ci rcuits 

From  the  discussion  in  section  3,  it  is  apparent  that  the  analysis 
of  lumped  nonlinear  circuits  reduces  to  the  repeated  analysis  of  the 
augmented  linear  circuit  excited  by  nonlinear  current  sources.  The 
terms  that  make  up  the  nonlinear  current  sources,  as  per  Table  3-1,  are 
factorable.  Thus,  the  "association  of  variables"  technique  can  be  ap¬ 
plied  to  obtain  the  transient  analysis  of  nonlinear  circuits.  We  now 
discuss  how  the  various  order  responses  can  be  obtained.  We  consider 
the  unit  step  input  case  here;  the  cases  for  other  factorable  inputs  can 
be  handled  similarly. 

FIRST  ORDER  RESPONSE; 

The  general  form  of  the  first-order  transfer  function  for  the  aug¬ 
mented  linear  network,  ignoring  multiple  poles*,  is: 


H. (s.)  =  £ 


P. 


VI 


i=i  Vpi 


,  M 

E 

k=Q 


R  sk 

Vi 


C4-33) 


where  P..,  p,.  and  R^  are  complex  constants.  Then  if  X(s^),  the  transform 
of  input  xlt),  is  factorable,  the  response  Y^Cs)  will  have  the  same  form 
as  eqn.  (4-33).  We  note  here  most  of  the  inputs  used  in  practice  are 
factorable;  that  is,  they  can  be  expanded  into  partial  expansion  form  as 

*The  case  of  multiple  poles  can  also  be  handled  by  the  "association 
of  variables"  technique.  In  our  discussion  here,  we  do  not  treat 
this  case. 
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in  eqn.  (4-33).  Specifically,  if  x(t>  is  a  step  function,  then  the  out¬ 


put,  Y^(s^),  will  have  the  form: 


A  N  A.  M-1  . 

Y.  (s.)=  —2.  +  £  — I — +  £  Q  sk 

1  1  S1  i=1  s1+pi  k=0  k  1 


(4-34) 


The  time-domain  response  can  be  easily  determined  from  eqn.  (4-34)  by 
taking  the  inverse  Laplace "transform. 

SECOND  ORDER  RESPONSE; 

The  equation  for  determining  the  second-order  transfer  function  is 
(see  section  3.3): 


H^s^s^  =  -Zfs^s^  I^s^s^  (4-35) 

where  is  the  vector  of  second-order  transfer  functions  for  the  multi¬ 
ple  node,  multiple  nonlinear  circuit,  Z(s^+s2)  is  the  inverse  of  the  re¬ 
duced  admittance  matrix,  and  ^(s^Sg)  is  the  vector  of  second-order 
current  sources.  For  the  four  types  of  nonlinear  elements  described  in 
section  3-3,  the  elements  of  the  vector  ^(s^s^  will  be  of  the  follow¬ 
ing  form,  depending  on  the  type  of  nonlinear  element: 

Nonlinear  Resistor: 

12<s1's2)  =  a2Hi(si>Hi(s2>  (4"36> 


Nonlinear  Dependent  Source: 


Ip(Sl/s2>  =  a^fs^H^Sj)  ♦  ag2H^(s^)H^(s2) 
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(4-37) 


+  \  •11ch^(s1>h^(s2>  +  H^(Si)h!J(s2)3 


Nonlinear  Capacitor; 

I2(S1'S2>  *  (s1+s2>a2H1(s1)H1<s2)  (4-38) 

Nonlinear  Inductor: 

I^(s1#s2>  =  a2H^(s1)H^(s2)/(s1+s2)  (4-39) 

where  I2  is  the  current  source  at  port  p,  a^  is  the  coefficient  of  the 
quadratic  term  of  the  nonlinearity,  and  H;J  is  the  first-order  transfer 
function  from  the  input  port  to  port  j.  Thus,  depending  on  the  number 
of  nonlinear  elements,  the  second-order  transfer  function  H*  can  be  ex¬ 
pressed  using  eqn.  (4-35)  as; 


H2(s1,s2 


K 

)  =£ 
k=1 


zlk<S1+S2>I2(s1/,S2> 


(4-40) 


where  z^ls^Sg)  is  the  l,k  element  of  Z(s,j+s2)  and  K  is  the  number  of 
second-order  current  sources.  - 

Observing  the  form  of  eqns.  (4-36)  through  (4-39),  we  recognize 
that  the  output  transfer  function  will  be  made  up  of  a  summation  of 
terms  of  the  following  form; 

H^s^Sj)  =  t(s1+s2)H1(s1)H1(s‘2)  (4-41) 

where  fl2  denotes  the  "partial"  transfer  function,  and  1  is  obtained 
after  the  (s^+s.^)  term  appearing  in  eqn.  (4-38)  and  (4-39)  has  been 
lumped  with  z(s^+s2).  Then  for  an  input  x(t)  having  transform  X(s),  we 
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have,  the  partial  output,  Y^Cs^s^),  given  by: 

=  fl2(s1,s2)X(s1>X(s2> 

=  i(s1+s2>H1(s1>X(s1>H1Cs2)X(s2>  (4-42) 

Recalling  that  H^(s^)XCs^)  is  the  linear  response  at  a  particular  port 
in  the  circuit,  we  can  re-write  eqn.  (4-4?)  as: 

?2(s1,s2)  =  J(s1+s2)Y1<s1)Y1(s2> 


!lSl*Sz)U  Vpi 


(4-44) 


Equation  (4-44)  has  the  same  form  as  eqn.  (4-30).  Therefore,  we  can  ap¬ 
ply  the  "association  of  variables"  technique  by  inspection  to  solve  for 
?2<s),  and  then  for  ?2<t).  By  taking  all  the  terms  in  eqn.  (4-40),  we 

can  obtain  y2(t),  the  second  order  time-domain  response  to  an  input 

k 

x(t).  It  should  be  noted  that  the  terms  2,  Q.  s  ,  which  were  present  in 

k  * 

eqn.  (4-34),  have  been  dropped  from  eqn.  (4-44).  This  is  because  func¬ 
tions  such  as  impulses,  doublets,  etc.  do  not  exist  when  squared. 


THIRD  ORDER  RESPONSE: 

The  equation  for  determining  the  third-order  transfer  function  vec¬ 
tor,  H3(s1,s2,s3),  is: 

-3<S1/S2'S3*  =  ■?(s1+82+S3>-3<S1/S2'S3)  (4-45) 

where  Z(*)  is  the  inverse  of  the  reduced  admittance  matrix  and  IT(*>  is 
the  third-order  current  source  vector.  The  elements  of  the  third-order 
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current  source  vector  depends  on  the  type  of  the  nonlinearity.  The  form 
of  these  elements  are  as  follows: 


Nonlinear  Resistor: 


lP(Sl,s2,s3>  =  2a2Hj(Sl)HP(s2/s3)  +  a^ls^H^s^H^)  (4-46) 


Nonlinear  Dependent  Source: 


lP<Sl,s2,s3>  =  2a2QH2(srs2)H!J<s3)  +  Za^H^s^H^s^Sj) 


+  a11(;Hj(s1)H^(s2,s3)  +  H2(srs2)H^(s3>: 


+  a30  .n  VV  +  aQ3  .n  »i<V 


(4-47) 


Nonlinear  Capacitor: 


I3(s.j,s2,s3)  =  (s1+s2+s3)C2a2H^(s1)H2(s2,s3) 


+  a3H^(s1)H^(s2)H^(s3)3 


(4-48) 


Nonlinear  Inductor: 


I^s^s^Sj)  =  C2a2HP(s1)HP(s2,s3) 
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'J  c  •,  -  : 

its.-.?*. 


■  ■  <  ■  r  7 


(4-49) 


+  a3H^(s1)HP(s2H^<s3)3/(s1+s2+s3) 

where  13<*)  is  the  current  source  at  port  p,  a^  is  the  coefficient  of 
the  nonlinearity  at  port  p,  and  and  H2  are  the  first-  and  second- 
order  transfer  functions  at  port  i.  Using  eqn.  (4-45),  the  third-order 
transfer  function  at  port  l  can  be  expressed  in  terms  of  the  individual 
third-order  current  sources  in  the  network  as  follows: 


H3<S1'VS3>  =  'u'VW's'VW 


(4-50) 


where  z^  is  the  l,k-th  element  of  matrix  Z  and  K  is  the  number  of 
third-order  current  sources. 

Observing  the  basic  form  of  eqns.  (4-46)  through  (4-49),  we  see  a 
striking  similarity.  There  are  basically  two  types  of  terms  making  up 
the  third-order  current  sources:  1)  due  to  the  squaring  operating,  giv¬ 
ing  rise  to  (s^H^s^Sj);  and  2)  due  to  the  cubing  operation  giving 

rise  to  n  H(s^).  The  transfer  function  at  port  l,  from  eqn.  (4-50), 
i=1  1 

is  a  summation,  a  typical  term  of  which  will  have  the  following  form: 
A3(s.j,s2/s3)  =  2(s1+2+s3)H<j(s1)H2(s2,s3) 


+  2(s1+s2+s3)H,|(s1)H1(s2)H1(s3)  (4-51) 

where  A3  is  the  "partial"  third-order  transfer  function,  and  2  is  ob¬ 
tained  after  the  (s^+s2+s3)  term  appearing  in  eqn.  (4-48)  and  (4-49)  has 
been  lumped  with  z(s^+s2+s3) . 

For  an  input  x(t)  having  a  transform  X(s),  the  partial  response, 
Y(s^,s2,s3),  Is  given  by: 
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=  fl3(s1/s2/s3)X<s1>X(s2)X(s3> 


=  C2(s1+s2+s3>H1 (s1)H2(s2/s3> 

3 

+  2(s1+s2+s3)H1(s1>H1(s2)H1(s3)D  n  X(sJ  (4-52) 

We  now  seek  to  obtain  ?3(s),  the  transform  of  93(t),  from  eqn.  (4-52). 
To  do  this,  we  look  at  each  term  in  eqn.  (4-52)  separately. 

Recalling*  that  H^s^s^  =  Eils.j+s^H.j  (s^H.j  (s2),  we  can  re-write 
the  first  term  of  eqn.  (4-52)  after  performing  the  symmetrization  opera¬ 
tion  as: 


2(s1+s2+s3)2(s2+s3)H1 (s1 )H1 Cs^H,, (s3)X(s1)X(s2>X(s3) 

=  ^4(s1+s2+s3)E2(s1+s2>+2(s1+s3>+2(s2+s3>3 
3 

n  H.(s . )X(s . )  (4-53) 

i=1  1  1  1 

Equation  (4-53)  is  in  a  form  suitable  for  applying  the  “association 
of  variables"  technique.  We  associate  two  variables  at  a  time  to  reduce 
?3(s^,s2,s3)  to  ?j(s).  For  example,  consider  the  following  association 
steps: 


♦see  eqns.  <4-40)  and  (4-41). 
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2(s1+s2»Sj)2(s^+s2)  [I  H1 (s..)X(s ..) 


=  2<s.j+S2+s3)  2<s1+s2>Y^  ^s1^1  (s^) 


Ss1's2 


=  2(s  +s,)Q(s  ) Y.  (s,) 
n  3  n  T  3 


aVS3 


=  R(s) 


(4-54) 


The  association  of  two  variables  is  identical  to  the  association  in¬ 
volved  in  determining  the  second  order  response.  Thus,  the  partial 
third-order  response  is  obtained  by  repeating  this  step  once. 

The  association  of  variables  involved  in  the  second  term  in  eqn. 
(4-51)  is  quite  straightforward.  Recalling  that  H^fs^Xfs.)  =  ^(s.)  we 
can  write  the  second  term  in  eqn.  (4-51)  as: 


2(s.+s?+s,)  n  Vs.)  =  2(s.+s-+s,)  n  (  -  L-  ■) 
1  c  5  i*1  1  1  1  * ■  *  i=1  j  si  pj 


(4-55) 


A  typical  term  in  eqn.  (4-55)  will  be  of  the  following  form: 


A,  A-  A_ 

P(s. ,s9,s,)  =  lls.+s.+s,)  — j--  .  tt 

123  123  Si+Pf  s2+p2  S3+P3 


(4-56) 


Associating  two  variables  at  a  tine,  the  final  form,  P(s),  is  given  by: 
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(4-57) 


PCs)  =  2(s) 


A1A2A3 

s+p1+p2+p3 


In  our  above  discussion  we  isolated  the  typical  terms  present  .in 
the  higher  order  responses  to  illustrate  the  application  of  the  associa¬ 
tion  of  variables  techniques  to  reduce  these  terms  into  a  form  suitable 
for  obtaining  the  transient  response.  By  summing  up  the  partial 
responses,  we  can  obtain  the  complete  response  for  the  network. 

The  association  of  variables  involved  in  determining  the  complete 
response  can  be  done  by  inspection,  provided  we  have  all  information  re¬ 
garding  the  poles,  and  the  associated  residues,  of  the  linearized  sys¬ 
tem.  For  adapting  this  scheme  on  the  computer,  an  algorithm  which  pro¬ 
vides  all  this  information  must  be  used.  A  semi-symbolic  analysis  C293 
of  the  linearized  system  must  therefore  be  part  of  the  algorithm  for 
adapting  Volterra  series  computer  aided  analysis.  The  main  problem  in 
implementing  this  scheme  on  the  computer  is  the  need  to  repeatedly  per¬ 
form  a  partial  fraction  expansion.  For  nonlinear  resistive  networks, 
this  task  is  not  so  difficult  if  we  make  use  of  some  of  the  basic  pro¬ 
perties  of  lumped  linear  systems.  We  will,  however,  defer  this  topic 
until  chapter  5. 

Finally,  it  is  noted  that  the  zero-input  time  domain  response  can 
be  obtained  by  using  the  ideas  presented  above.  For  a  circuit  with  mul¬ 
tiple  dynamic  elements,  this  reduces  to  the  multiple-input,  multiple- 
output  problem,  where  each  input  corresponds  to  the  initial  condition 
across  the  energy  storage  element. 


Ill 


CHAPTER  5 


COMPUTER-AIDED  ANALYSIS  USING  VOLTERRA  SERIES 

5-1 .  Introduction 

The  adapting  of  Volterra  series  method  in  a  general  simulation  pro¬ 
gram  has  been  regarded  as  difficult  by  various  authors  116,29,43].  As 
such,  virtually  no  effort  has  been  spent  on  investigating  the  computa¬ 
tional  aspect  of  this  method  and  using  it  for  the  spectrum  and  time- 
domain  analysis  of  general  nonlinear  circuits  with  polynomial  type  non- 
linearities.  Previous  works  [7,36]  have  endeavored  to  check  the  validi¬ 
ty  of  this  approach  by  applying  it  to  specific  circuit  problems  using  a 
computer,  but  have  never  implemented  the  approach  in  a  general  simula¬ 
tion  program. 

The  only  major  effort  in  using  Volterra  series  for  general  non¬ 
linear  circuit  analysis  has  been  the  development  of  the  program  NCAP 
£10,44],  A  cursory  review  of  this  program  reveals  the  inherent  ineffi¬ 
ciency  in  the  computational  approach  with  regards  to  storage  and  types 
of  algorithms  used.  This  inefficiency  notwithstanding,  there  are  severe 
limitation  regarding  the  usefulness  of  the  program:  first,  the  program 
merely  computes  the  numerical  values  of  the  nonlinear  transfer  function 
at  the  various  program-prescribed  combinations  of  the  input  frequencies, 
and  does  not  compute  all  the  transfer  function  values  which  are  required 
to  compute  the  complete  output  spectrum.  Thus,  NCAP  does  not  yield  the 
entire  output  spectrum  information.  Second,  to  compute  up  to  an  n-th 
order  transfer  function,  the  user  must  specify  n  input  frequencies, 
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which  are  assumed  to  be  a  sum  of  exponentials  and  not  real  sinusoids. 
The  program,  therefore,  is  severely  limited  in  its  usefulness  from  the 
point  of  view  of  a  user,  who  may  only  be  interested  in  obtaining  the 
output  spectrum  -  say,  for  example,  up  to  the  third  order  response  to 
two  sinusoidal  inputs  -  and  has  little  use  for  the  numerical  values  of 
the  transfer  functions  at  the  program  prescribed  frequencies. 

The  zero-state  time-domain  analysis  of  nonlinear  circuits  using 
Volterra  series  approach  on  a  computer  has  not  been  attempted  before  - 
not  even  for  a  specific  nonlinear  circuit  problem. 

In  this  section  we  look  at  the  computational  aspect  of  the  Volterra 
series  method  for  general  simulation  purposes  and  then  present  the  basic 
algorithms  for  adapting  this  method  for  1)  spectrum  and  distortion 
analysis,  and  2)  for  zero-state  time-domain  analysis  to  a  step  input. 
A  digital  computer  program,  PRANC  (Program  for  Analysing  Nonlinear 
Circuits),  which  makes  use  of  the  algorithms,  has  been  written  and  im¬ 
plemented  on  the  CDC  6500  computer.  A  detailed  description  and  program 
listing  is  contained  in  a  separate  technical  report  C41].  Some  examples 
from  the  use  of  this  program  are  given  in  section  6. 

In  section  5-2,  we  present  a  brief  overview  of  symbolic  analysis  in 
linear  circuits,  and  then  describe  the  reason  why  a  symbolic  approach  is 
particularly  useful  in  adapting  Volterra  series  for  general  simulation. 
Section  5-3  deals  with  the  implementation  of  the  symbolic  approach,  and 
also  contrasts  the  computational  effort  between  a  numerical  approach  and 
the  particular  symbolic  approach  used  here.  The  algorithm  for  obtaining 
the  complete  output  spectrum  and  the  various  distortion  indices  is 
described  in  section  5-4;  the  algorithm  for  zero-state  time-domain 
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analysis  to  a  step  input  is  presented  in  section  5-5.  A  description  of 
the  computer  implementation  of  these  algorithms  is  given  in  section  5-6. 


5-2.  Why  a  Symbolic  Analysis  Approach. 

The  symbolic  analysis  of  circuits  invoLves  the  computation  of  the 
a.  and  b.  for  network  functions  in  the  form 

i  i 


F<s>=Me.  V‘i  — 


D(s) 


Eb.  s1 


(5-1) 


when  all  circuit  elements  are  known.  The  more  general  form 


F(s;x.,x_,...,x  )  • 
i  c  n 


N(s;x^. . -/X^) 

D (s;x. , . . ./X  ) 
i  n 


(5-2) 


applies  when  some  elements  of  the  circuit  x^  are  kept  as  symbols.  The 
advantages  of  symbolic  analysis  have  been  recognized  previously  C29,45J. 
One  particular  advantage,  and  the  one  which  is  relevant  to  our  problem 
here,  is  that  the  numerical  evaluation  of  a  function  at  discrete  points 
is  much  easier  and  faster  once  the  symbolic  function  is  obtained  than 
working  repeatedly  with  a  circuit  analysis  program.  With  this  brief 
overview  of  symbolic  analysis,  we  now  proceed  to  answer  the  question: 
Why  use  a  symbolic  analysis  approach  for  adapting  the  Volterra  series 
method  for  general  circuit  analysis? 

As  pointed  out  in  the  previous  sections,  a  nonlinear  circuit  is 
completely  characterized  by  its  Volterra  kernels,  or  their  transforms  - 
the  generalized  transfer  functions.  These  transfer  functions  are  then 
directly  related  to  the  various  order  sinusoidal  steady-state  responses, 
as  described  in  Chapter  4.  The  n-th  order  transfer  function  is  deter- 
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mined  from  the  following  equation  (see  Chapter  3): 


H  (s.,,...,s  )  =  CY< 
-nr  n 


s..)3-1  I 


Is* 


i=1 


-n  1 


'V 


(5-3) 


n 

where  Y(  £  s.)  is  the  reduced  node  admittance  matrix  evaluated  as  s. 

~  i=1  1  1 

+  s->  +  ...  +  s  ,  and  I  is  the  n-th  order  current  source  vector  due  to 

2  n  -n 

the  nonlinear  elements.  To  compute  the  output  spectrum,  we  evaluate  Hn 
at  various  and  many  frequency  combinations.  From  eqn.  (5-3)  it  should 
be  clear  that  such  an  evaluation  will  entail  the  inversion  of  the  re¬ 
duced  node  admittance  matrix  at  each  of  these  frequency  combinations. 
Using  combinational  analysis,  it  has  been  shown  [403  that  for  an  input 
consisting  of  M  sine  waves,  the  number  of  inversions  involved  in  an  n-th 

order  response,  given  by  N  ,  is: 

n,m 


N 


n,m 


(5-4) 


Thus,  for  a  5-tone  input  and  up  to  a  third  order  analysis,  the  number  of 
inversions  is  approximately  285.  For  higher  order  responses,  this  num¬ 
ber  grows  up  very  rapidly. 

There  are  two  basic  approaches  available  to  handle  this  inversion 
process:  1.  Numerical  approach,  or  2.  Symbolic  approach.  The  advan¬ 
tage  of  evaluating  symbolic  transfer  functions  mentioned  earlier  makes 
the  symbolic  approach  more  attractive.  How  much  advantage  is  gained  in 
using  a  symbolic  analysis  depends  on  how  much  computational  effort  is 
expended  in  obtaining  the  symbolic  inverse  of  the  reduced  node  admit¬ 
tance  matrix;  thus,  we  need  an  efficient  scheme  for  obtaining  the  sym¬ 
bolic  inverse.  The  determination  of  the  symbolic  inverse  will  be  the 
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subject  of  section  5-3. 

Another  reason  for  using  a  symbolic  analysis  is  concerned  with  the 
determination  of  the  zero-state  time-domain  response.  The  basic  step 
involved  in  the  "association  of  variables"  technique,  which  was  present¬ 
ed  in  section  4-3,  is  the  prescribed  combining  of  the  various  poles,  and 
their  associated  residues,  of  the  augmented  linear  network.  A  pole  lo¬ 
cation  and  its  residue  can  be  readily  determined  once  the  transfer  func¬ 
tion,  as  per  eqn.  (5-1),  is  known.  The  entries  of  the  symbolic  inverse 
of  the  reduced  node  admittance  matrix,  which  form  the  higher  order 
transfer  functions,  will  contain  all  the  information  needed  to  perform 
the  prescribed  combining  of  the  various  poles  and  residues  to  determine 
the  zero-state  response. 

The  reasons  presented  above  stem  from  looking  at  the  computational 
aspect  of  adapting  Volterra  series  for  computer-aided  analysis.  There 
are  other  advantages  gained  from  using  a  symbolic  analysis.  An  impor¬ 
tant  one  is  that  the  generalized  transfer  functions  can  be  obtained  as 
functions  of  s.  once  the  inverse  of  the  reduced  node  admittance  matrix 
is  obtained  as  a  symbolic  function  of  s.  This  can  be  seen  from  examin¬ 
ing  eqn.  (5-3).  The  formation  of  the  n-th  order  current  source  vector 
is  a  bootstrapping  operation,  as  was  pointed  out  in  Chapter  3.  That  is, 
an  n-th  order  source  is  formed  from  transfer  functions  of  order  less 
than  n.  The  first-order  transfer  function  vector  is  determined  from  a 
column*  of  the  symbolic  inverse  of  the  reduced  node  admittance  matrix. 
The  second  order  current  sources,  which  depend  on  the  elements  of  the 
first  order  transfer  function  vector,  are  therefore  formed  from  this 

*This  is  assuming  a  single  input  circuit. 


column  of  CY(s)3  \  The  second-order  transfer  function  vector  is  ob¬ 
tained  by  pre-multiplying  the  second-order  current  source  vector  by 
-1 

CYls^+s^)]  ,  according  to  which  the  second-order  transfer  function  vec¬ 
tor  eventually  depends  on  the  entries  of  inverse  of  the  node  admittance 
matrix  evaluated  at  (s^+s^).  The  third-  and  higher-order  transfer  func¬ 
tions  have  a  similar  dependence.  Thus,  if  the  inverse  of  the  reduced 
node  admittance  matrix  is  obtained  in  symbolic  form,  with  s  retained  as 
a  symbol,  then  a  functional  description  of  the  nonlinear  functions  can 
also  be  obtained.  A  concomitant  advantage  of  this  functional  descrip¬ 
tion  is  that  theorems  from  multi-dimensional  theory  L5D  (such  as  initial 
value,  final  value,  etc.)  can  then  be  applied  to  these  transfer  func¬ 
tions  to  gain  more  insight  into  the  workings  of  the  circuit. 

In  Chapter  2  we  developed  recursive  relationships  to  estimate  the 
error  incurred  in  the  truncation  of  the  series  solution.  This  error  was 
directly  related  to  the  l^  norm  of  the  linear  kernel  function,  which,  in 
turn,  is  related  to  the  poles  and  residues  of  the  linearized  system. 
Thus,  we  can  get  an  estimate  of  the  accuracy  of  our  solution  through  the 
pole-residue  information  provided  to  us  by  the  symbolic  analysis. 


5-3.  Symbolic  Analysis  Method 

Symbolic  circuit  analysis  by  digital  computer  has  been  of  consider¬ 
able  interest  in  the  past  decade.  Many  algorithms  and  methods  have  been 
derived  to  obtain  symbolic  transfer  functions  of  linear  circuits  C29]. 
Most  of  these  methods  use  tree  enumeration  C463,  signal-flow  graphs 
Z291,  or  purely  numerical  methods  C473  to  obtain  symbolic  transfer  func¬ 
tion  between  the  input  and  the  output.  These  approaches  are  basically 
useful  for  single-input,  single-output  systems.  The  inversion  of  the 


reduced  node  admittance  matrix  to  obtain  the  open-circuit  impedance  ma¬ 
trix,  which  is  the  problem  we  are  dealing  with,  is  basically  a  multi¬ 
input,  multi-output  problem.  The  methods  mentioned  above  can  be  adapted 
to  solving  the  problem  at  hand;  however,  the  generation  of  multiple  sym¬ 
bolic  functions  using  these  approaches  many  not  be  satisfactory  because 
of  excessive  computer  time  requirements.  Some  other  approach  is  defin¬ 
itely  warranted. 

Published  methods  L25-273  for  inverting  the  nodal  admittance  matrix 
when  the  elements  are  rational  functions  of  the  Laplace  transform  vari¬ 
able  s  use  pivotal  techniques.  It  may  appear  that,  since  it  is  easy  to 
program  a  computer  to  perform  polynomial  arithmetic,  these  pivotal- 
techniques  are  a  natural  way  to  approach  the  symbolic  inversion  problem. 
Results  from  the  use  of  such  a  technique  have  proved  to  be  disappoint¬ 
ing,  mainly  due  to  the  following  reasons: 

(a)  The  process  of  inversion  transforms  the  nodal  admittance  matrix, 

which  contains  terms  of  the  form  as  +  ~  +  c,  into  a  matrix  in  which 

s 

every  element  is  a  rational  function  of  s.  The  pivotal  technique  pro¬ 
duces  the  inverse  matrix  where  common  factors  appear  between  numerator 
and  denominator,  and  unless  some  mechanism  is  built  into  the  process 
whereby  these  common  factors  are  recognized  and  removed,  the  elements 
produced  will  have  polynomials  of  excessively  high  order. 

(b)  When  the  circuit  complexity  is  high,  the  evaluation  of  the  symbolic 
function  at  high  frequency  values  can  give  rise  to  numerical  problems. 

g 

For  example,  a  circuit  with  8  poles  will  have  an  s  term  in  the  charac¬ 
teristic  polynomial.  When  evaluated  at  10  Mrad/sec,  this  term  produces 
a  number  equal  to  10^.  Of  course,  this  problem  can  be  alleviated  by 
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obtaining  a  partial  fraction  expansion  (PFE)  form  for  the  transfer  func¬ 
tions.  But  this  again  entails  additional  computations  -  not  to  mention 
the  numerical  instability  problems  involved  in  root  finding. 

(c)  It  has  also  been  found  that  pivotal  techniques  become  numerically 
unstable  for  higher  order  circuits. 

We  therefore  seek  another  alternative  for  obtaining  the  symbolic 
open  circuit  impedance  matrix. 

An  approach  based  on  the  state  variable  formulation  can  be  used  to 
achieve  this  goal.  Specifically,  consider  the  general  p-port  augmented 
linear  circuit  of  Fig.  5-1 (a).  We  wish  to  solve  for  the  transfer  im¬ 
pedances,  2^(s),  i,j  =  1,2,...,p,  from  the  j-th  port  to  the  i-th  port. 
Knowing  these  transfer  impedances,  we  can  write  for  the  p-port: 

V(s)  =  Z(s)I(s)  *  CY(s)3_1  I(s>  (5-5) 

—  — 


where  V(s)  =  CV„ (s)  V,<s>  ...  V  (s)D 

-  1  2  p 


(5-6) 


Z(s)  =  Cz. . (s)3 
~  i) 


(5-7) 


and 


Hs)  =  CI1  (s)  I2(s) 


I  (s)D 
P 


(5-8) 


Note  that  the  vector  V(s)  contains  entries  which  are  the  output  voltages 
and  voltages  that  control  the  nonlinear  element  characteristics  in  the 
nonlinear  circuit. 

To  obtain  Z(s)  symbolically,  we  write  for  the  network  of  Fig. 
5-1 (b),  the  following  state  equations: 
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l  =  A* +  Si 


(5-9) 


v  =  Cx  (5-10) 

where  £  is  the  vector  of  state  variables,  and  v_  and  j_  are  vectors  whose 
transforms  appear  in  eqns.  (5-6)  and  (5-8),  respectively.  Taking  the 
Laplace  transform  of  eqn.  (5-9)  and  (5-10),  and  solving  for  V(s),  we 
get: 

V(s)  =  CC(sI-A)"1  B  *  03  _I(s)  (5-11) 

and,  therefore,  we  get  Z(s>  to  be 

Z(s)  =  CC(sI-A)"1  B  +  03  (5-12) 

which  is  identically  the  inverse  of  the  reduced  node  admittance  matrix. 
The  matrix  (sI-A)  can  be  inverted  by  applying  the  similarity 

•V  <V 

transformation  as  follows: 

A  =  M  a  M_1  or  a  =  M_1  A  M 

M_1  (sI-A)M  =  si  -  H_1AM  =  si  -  A 

or  (Sl-A)'1  =  M(sl-A)'1  M'1  (5-13) 

-1  -1 

where  the  inverse  of  (si  -  a)  is  simply  ,diag  ■C(s-x^)  ,(s-X2J  r  •••* 
where  x^  are  the  eigenvalues*  of  the  A  matrix  and  H  is  the  modal  matrix. 
Substituting  eqn.  (5-13)  Into  eqn.  (5-12),  we  get, 

*Here  we  assume  distinct  eigenvalues;  the  repeated  eigenvalues  can 
be  handled  similarly. 
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(5-14) 


Z(s)  =  CCM(sI-a)"1  M_1  B  +  DD 

=  CC(sI-a)'1  6  +  D3 

where  £  £  CM  and  S  |  M  ^B.  Equation  (5-14)  yields  the  entries  of  Z(s) 
in  partial  fraction  expansion  form,  which,  as  mentioned  previously,  is  a 
more  desirable  form  from  a  computational  standpoint.  All  information 

A  A 

regarding  Z(s)  is  contained  in  the  matrices  C,  B,  D  and  a  vector  con¬ 
taining  the  eigenvalues.  An  algorithm  for  implementing  this  approach  is 
given  in  Fig.  5-2.  It  should  be  noted  that  the  approach  used  here  is 
completely  numerical  and  does  not  involve  any  coding  and  decoding  of 
symbols. 

Now  that  an  algorithm  for  obtaining  the  symbolic  Z(s)  is  defined, 
we  can  make  a  comparison  of  the  computational  effort  involved  between 
using  a  symbolic  inverse  and  the  numerical  inverse  of  the  node  admit¬ 
tance  matrix  at  each  frequency  point. 

The  computational  trade-off  between  the  symbolic  approach  and  a  nu¬ 
merical  approach  for  matrix  inversion  is  very  problem  dependent.  While 
a  clear-cut  winner  cannot  be  established,  a  tentative  answer  can  be  ob¬ 
tained  by  noting  the  operations  count,  defined  in  terms  of  multiplica¬ 
tions  and  additions,  involved  in  the  two  schemes. 

In  the  case  of  the  numerical  approach,  the  number  of  independent 
nodes,  n,  and  the  number  of  branches,  b,  are  the  most  important  quanti¬ 
ties  for  determining  the  computational  effort  along  with  the  number  of 
frequency  points  at  which  the  output  is  desired.  Assuming  that  no 

sparse  matrix  techniques  are  used,  the  numerical  inversion  of  an  (nxn) 

* 

matrix  requires  Q(n  /3>  units  of  work,  where  0(  )  *  "order  of",  and  1 
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Fig.  5-2.  Algorithm  for  inverting  Y(s>  symbolically 


unit  of  work  =  one  addition  and  one  multiplication.  For  k  frequency 
points,  the  work  becomes  0(kn^/3).  This  does  not  involve  book-keeping 
and  other  pre-  and  post-processing  steps  such  as  pivoting  and  iterative 
refinement,  which  are  usually  necessary  to  insure  reliability  and 
robustness  of  the  algorithm. 

In  the  case  of  symbolic  inversion  using  our  approach,  the  important 
parameters  in  the  computational  effort  are  the  dynamic  degrees  of  free¬ 
dom,  d,  and  the  number  of  ports,  p,  where  voltages  and  currents  are  in¬ 
jected  or  measured.  Using  the  QR  algorithm  £29,483  for  computing  the 
eigenvalues  of  the  A  matrix,  the  operation  count  is  0(8d^).  The  total 
work  required  for  obtaining  the  inverse  at  k  frequency  points  is  there¬ 
fore  0(8d^  +  kdp^).  The  number,  p,  depends  on  the  number  of  nonlineari¬ 
ties  in  the  circuit,  and  is  usually  small.  Also,  if  the  network  com¬ 
plexity  is  less  than  the  number  of  nodes,  the  symbolic  approach  would, 
in  general,  require  less  computational  effort.  As  far  as  accuracy  is 
concerned,  both  the  QR  algorithm  and  the  Crout's  algorithm  with  pivoting 
and  iterative  refinement  yield  accurate  results. 

The  efficiency  of  the  symbolic  method  rests  heavily  upon  the  avai¬ 
lability  on  an  efficient  process  for  forming  the  state  equations.  The 
hybrid  analysis  method  £28,293,  which  essentially  reduces  to  the 
analysis  of  a  resistive  network,  is  well-suited  for  our  purposes  here. 
We  shall  discuss  this  topic  in  section  5-6. 

5-4 •  Spectrum  and  Distortion  Analysis  Algorithm 

The  output  spectrum  and  distortion  indices  for  a  nonlinear  circuit 
with  polynomial  type  nonlinearities  can  be  computed  on  the  basis  of  the 
material  of  Chapters  3  and  4.  A  flow-chart  of  the  basic  algorithm  for 
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such  a  computation  is  given  in  Fig.  5-3.  We  describe  the  steps  involved 
in  the  following  paragraphs: 

Step  1 :  For  the  given  nonlinear  circuit,  determine  the  dc  operating 
point.  Expand  each  nonlinear  function  into  a  Taylor  series  about. the 
operating  point  to  get  a  polynomial  representation  for  the  nonlinear 
element  in  terms  of  the  incremental  quantities.  Thus,  for  example,  a 
forward-biased  diode  having  the  "global"  V-I  representation 

I  =  IsCexp(qV/nkT)  -  13  (5-15) 

can  be  expanded  into  a  Taylor  series  to  yield  the  following  incremental 
v-i  representation: 

1  =  xo  ISkT  v  +  7T  (7JkT)2  *  *  3T  (7*t>3  *  +  —  (5‘16> 

where  Ig  is  the  dc  operating  current. 

Step  2:  Lump  the  linear  part  of  the  nonlinear  elements  with  the  existing 
linear  network  to  form  the  augmented  linear  network.  Extract  as  ports 
the  nonlinear  element  branches  and  the  branches  that  control  the  non¬ 
linear  clement  characteristics  (dependent  nonlinear  element  case),  along 
with  the  output  and  source  branches,  from  the  augmented  linear  network. 
Let  V  =  CV.  V,  ...  V  3  and  I  *  Cl.  I-  ...  I  3  denote  the  vector  of  vol- 
tages  and  currents  for  these  ports,  respectively. 

Step  3:  Using  a  symbolic  analysis  algorithm  (see  Fig.  5.2),  obtain  the 
entries  of  the  Z  matrix  as  a  function  of  s,  where 

V(s>  *  Z(s>  l<s>  (5-17) 

For  each  of  the  input  sources,  and  their  associated  frequency  tones. 
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Figure  5-3.  Algorithm  for  Spectrum  and  Distortion  Analysis. 
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compute  the  first-order  output  voltages  at  each  of  the  extracted  ports 
by  using  the  appropriate  entries  of  the  Z  matrix.  This  step  amounts  to 
letting  s  =  j«>.  in  z^.Cs),  the  entries  of  Z^s). 

Step  4:  The  second-order  output  spectrum  is  evaluated  using  the  follow¬ 
ing  relationship: 

V2ls^,S2>  =  Z(s^+S2)  I^s,.,^)  (5-18) 

The  vector  J^(s,j,S2)  is  the  second-order  current  source  vector,  which  is 
formed  by  using  the  coefficients  associated  with  the  quadratic  term  of 
the  nonlinear  element  and  the  first-order  output  at  the  controlling 
port(s)  of  the  nonlinearity.  The  latter  information  was  obtained  in 
step  3.  The  given  input  tones  are  taken  two  at  a  time  in  eqn.  (5-18), 
along  with  the  information  derived  in  section  4-2,  to  evaluate  the  out¬ 
put  voltages  at  each  of  the  p-ports. 

The  third-order  output  spectrum  is  obtained  in  exactly  the  same 
manner.  The  first-  and  second-order  outputs  are  used  to  form  the 
third-order  current  source  at  each  combination  frequency,  which  is  then 
pre-mult ip  lied  by  evaluating  Z(s)  at  the  combination  frequency. 

Step  5:  Perform  a  histogram  analysis  of  all  frequency  points  and  combine 
the  responses  at  points  which  are  repeated.  The  distortion  indices  are 
computed  using: 

|V0(2«.)| 

2  =  TV^TT 


HD 


(5-19) 


(5-20) 


|V  (3<0| 

hd3  =  iv0ror 


where  HO^  and  HD^  denote  the  second  and  third  order  harmonic  distortion 
indices. 


5-5 .  Time-Domain  AnaLysi s  Algorithm. 

With  Z(s)  in  eqn.  (5-17)  obtained  in  symbolic  form,  the  zero-state 
time-domain  response  to  Laplace  transformable  inputs  can  be  determined 
by  extending  the  algorithm  discussed  in  section  5-4.  This  extension  is 
carried  on  the  basis  of  the  material  of  section  4-3.  Before  outlining 
the  new  algorithm  here,  we  derive  relationships  for  obtaining  the  com¬ 
plete  second-  and  third-order  responses  to  a  step  input.  The  case  of 
nonlinear  capacitor  and  inductor  is  not  treated  here.  As  mentioned  in 
C29],  any  nonlinear  network  can  be  equivalently  represented  by  a  non¬ 
linear  resistive  network;  we  therefore  treat  this  case  hert. 

Consider  a  network  with  multiple  nonlinear  resistive  elements,  as 
shown  in  Fig.  5-4.  As  per  steps  2  and  3  of  section  5-4,  we  have  the  re¬ 
lationship  (5-17)  for  the  augmented  linear  network.  Then,  for  a  lumped 
linearized  system,  the  transfer  function  H^j  (s),  i=1,2,...,p,  is  of  the 
following  form*: 


...  V.(s) 

V  (s)  =  v’  (s)  s  2i1(s>  = 
in 


N 

£ 

j=i 


(il ) 


1 


S+Pj 


M 

£ 

k=0 


(5-21 ) 


Then,  for  a  step  input,  the  first-order  output  at  each  of  these  ports 
is: 

*Without  any  lose  of  generality,  we  have  assumed  that  the  source 
port  is  the  port  number  1. 
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(5-22) 


(5-25) 


N 

<  E 


k=0 


a(l) 

3_) 

s2+p 


Y2^(s.j,s2)  can  be  reduced  to  Y2^(s),  the  transform  of  y2^(t),  by  us¬ 
ing  the  association  of  variables  C53.  This  yields  the  following: 


Y^Cs)  = 


L 

-  E 

1=1 


<l)r, 
a2  C( 


N 

E 

d=1 


(il) 


S+Pj 


-)  < 


N 

E 

n=0 


N 

E 

m=0 


(D  .(D 


n 


m 


S+P_+P_ 
n  m 


-) 


(5-28) 


Since  the  input  is  a  step  function,  eqn.  (5-28)  can  be  simplified  to*: 


v(i), 

Y2  (s) 


B(i)  r(i) 

£  !i_,  £  e* 

k  s+£k  j  (s+Pj )* 


(5-29) 


where  Bk  and  are  related  to  r^'s  and  A 
linear  system  poles,  p^ .  The  second-order  time  response  can  be  easily 
obtained  from  eqn.  (5-29). 

The  third-order  time  domain  response  can  be  obtained  similarly.  Me 
do  not  present  it  here  in  order  to  conserve  space.  Suffice  it  to  say, 
the  third-order  response  can  be  obtained  from  the  values  of  the  poles, 
and  their  associated  residues,  of  the  linearized  system  without  perform¬ 
ing  any  numerical  integration  or  iteration.  We  now  outline  the  algo¬ 
rithm  for  obtaining  the  zero-state  time-domain  response  of  a  network 
with  resistive  nonlinearities  (see  Fig.  5-4): 

Step  1:  From  the  polynomial  description  of  the  resistive  nonlinearities, 
form  the  augmented  linear  network.  Extract  as  ports  the  nonlinear  ele¬ 
ment  branches  along  with  the  output  and  the  source  branches.  Let 

*Here  we  have  assumed  that  the  linearized  system  had  distinct  poles. 
The  case  of  multiple  poles  can  be  handled  similarly,  but  it  compli¬ 
cates  the  algorithm  for  computer  implementation. 


k's,  and  ^  is  related  to  the 
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Fig.  5-4.  Algorithm  for  computing  time-domain  response  using 

Volterra  series. 
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V  =  CV.  V0  ...  V  3  and  I  =  Cl.  I,  ...  I  3  denote  the  vector  of  voltages 

—  \  d  p  —  I  £  P 

and  currents  for  these  ports,  respectively. 

Step  2;  Using  a  symbolic  analysis  algorithm  (see  Fig.  5-2),  obtain  the 
entries  of  the  Z  matrix  in  partial  fractial  expansion  form,  where 


V(S)  =  Z (s)  I(S) 
— 


(5-30) 


Step  3;  From  the  pole-residue  information  obtained  in  step  2,  obtain  the 
first-order  output  response  to  the  Laplace  transformable  input.  This 
amounts  to  evaluating  the  sum  of  exponentials  at  each  time  step. 

Step  4;  Next  obtain  the  second-  and  higher-order  responses  by  appropri¬ 
ately  associating  the  poles  and  residues  of  the  linearized  system. 
Thus,  for  example,  in  the  second-order  case,  this  amounts  to  determining 
Bk,  Cfe,  and  in  equation  (5-29).  Once  this  information  is  obtained, 
the  output  is  determined  by  evaluating  terms  of  the  form  Bk  expl^t)  and 
Ck  t  exp  (p^t)  at  each  discrete  time  point. 


5-6.  Program  PRANC . 

The  Program  for  Analysing  Nonlinear  £ircuits,  known  as  PRANC,  is  a 
digital  computer  program,  written  in  FORTRAN  IV,  that  computes  up  to  the 
third-order  complete  output  spectrum  of  a  nonlinear  circuit  with  polyno¬ 
mial  nonlinearities  driven  by  up  to  two  multi-frequency  inputs.*  In  the 
process  it  computes  the  Volterra  transfer  functions  at  each  of  the  fre¬ 
quency  combinations  involved.  PRANC  also  computes  the  zero-state  time- 
domain  response  of  circuits  with  only  resistive  nonlinear  elements  and  a 
step  input. 

*Thus,  mixer-type  circuits  can  be  analyzed  using  PRANC. 
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As  mentioned  previously/  the  solution  of  the  nonlinear  circuit 
problem  reduces  to  the  repeated  solution  of  the  linear  circuit.  To  ef¬ 
ficiently  handle  this  basic  problem/  PRANC  uses  a  semi -symbolic  approach 
L291  for  analysing  the  augmented  linear  circuit.  Specifically/  the  in¬ 
verse  of  the  reduced  node  admittance  matrix  is  obtained  in  terms  of  the 
symbol  s  using  the  state  equation  formulation  as  described  above. 

The  state  equations  for  the  linear  circuit  are  formulated  via  the 
Hybrid  analysis  method  C28/293.  If  T  denotes  port  branches  in  the  tree 
1291  and  C  denotes  port  branches  in  the  co-tree  of  a  linear  circuit/ 
then  the  Hybrid  analysis  yields  the  following  relationship: 


H11  H12  H13  H1 

H21  H22  H23  H2 


H3 

H, 


H2 

H. 


H2 

H, 


V 

V 

c 

VT 

i 

c 

=  0 


(5-31) 


■■  "  V  "  ■  '  " 

H  z 

By  suitably  forcing  the  various  ports  in  the  linear  circuit  into  the 
tree  and  co-tree  branches/  PRANC  uses  the  above  formulation  for  setting 
up  the  state  equations.  All  capacitor  branches  are  extracted  as  ports 
which  necessarily  become  part  of  the  tree  and  all  inductors/  nonlinear 
element  branches  (which  are  assumed  to  be  voltage  controlled)/  and  input 
and  output  branches/  are  extracted  as  ports  which  are  forced  as  part  of 
the  co-tree.  The  matrix  H  is  obtained  in  a  form  where  111  “  i  (I  being 
the  identity  matrix)/  ^  =  0,  H 22  »  I.  This  yields  the  capacitor 

currents  and  the  inductor  and  nonlinear  element  branch  voltages  in  terms 
of  known  variables.  Thus/  the  A/  B/  C,  and  D  matrices  in  the  state  and 
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output  equations  (see  eqns.  5-9  through  5-12)  are  obtained  from  the  sub¬ 
matrices  of  H.  The  formulation  of  eqn.  (5-31)  is  quite  fast,  since  it 
only  involves  the  analysis  of  a  resistive  network. 

It  is  noted  that  the  matrix  H  may  not  exist  in  idealized  circuits. 
However,  for  most  practical  circuit  this  matrix  is  almost  certain  to  ex¬ 
ist  [293.  It  should  also  be  noted  that  the  above  formulation  of  state 
equations  tacitly  assumes  that  no  degenerate  cutsets  (all  inductor- 
current  source  cutset)  or  degenerate  loops  (all  capacitor-voltage  source 
loop)  are  present  in  the  linearized  circuit.  These  restrictions  are  not 
very  severe,  especially  when  the  realistic  lossy  models  of  circuit  com¬ 
ponents  are  taken  into  account. 

The  next  step  in  the  PRANC  algorithm  is  to  determine  the  eigen¬ 
values  and  the  eigenvectors  of  the  A  matrix.  For  this  purpose,  the  dou¬ 
ble  QR  algorithm  [483  for  obtaining  the  eigenvalues  is  employed.  The 
basic  steps,  such  as  matrix  balancing,  reduction  to  Hessenberg  form, 
shift  of  origin,  are  employed  in  this  algorithm  to  make  it  efficient  and 
reliable.  The  eigenvectors  are  also  obtained  in  the  process. 

All  Information  about  the  inverse  of  the  reduced  node  admittance 

A  A 

matrix  is  stored  as  three  matrices  and  a  vector.  The  matrices  are  8,  C, 
and  D  (see  eqns.  5-14),  and  the  vector  contains  the  eigenvalues.  It  is 
noted  that  the  solution  of  eigenvectors  for  repeated  eigenvalues  can  be 
a  numerical  unstable  process  [493.  Thus,  the  programs  outputs  a  diag¬ 
nostic  message  when  such  a  case  occurs. 

The  first-order  voltage  response  at  the  prescribed  ports  is  now 
computed  from  the  entries  of  the  open-circuit  impedance  matrix.  These 
ports  Include:  source  port,  output  ports,  nonlinear  element  ports,  and 
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ports  which  control  the  nonlinear  element  characteristics.  The  response 
is  calculated  for  each  user  prescribed  frequency,  and  stored  as  a  two- 
dimensional  array:  port  number  vs.  the  frequency  number. 

The  second-order  voltage  response  is  computed  at  each  distinct  com¬ 
bination  of  the  input  tones  taken  two  at  a  time.  The  ports  of  interest 
are  the  same  as  that  for  the  first-order  response.  The  second-order 
current  source  vector,  at  a  particular  frequency  combination,  is  formed 
by  considering  the  nonlinear  element  type  and  the  voltage(s)  controlling 
it,  which  is  determined  from  the  first  order  response  array.  This  vec¬ 
tor  is  pre-mu Iti plied  by  the  open-circuit  impedance  matrix  evaluated  at 
the  combination  frequency  to  obtain  the  second-order  transfer  function 
vector  at  that  frequency.  The  response  voltage  at  this  frequency  is 
then  determined  from  the  transfer  function  value.  The  second-order 
transfer  function  values  are  again  stored  as  a  two-dimensional  array: 
port  number  and  the  particular  frequency  combination. 

The  third-order  response  is  determined  similarly.  The  third-order 
current  source  vector  is  formed  by  properly  picking  out  the  values  of 
the  first-  and  second-order  transfer  functions.  The  indexing  of  the  ar¬ 
rays  is  of  critical  importance  to  the  efficient  implementation  of  this 
scheme. 

When  the  time-domain  response  of  a  circuit  with  only  nonlinear 
resistive  elements  is  desired,  the  procedure  is  identical  up  until  the 
formation  of  the  symbolic  open-circuit  impedance  matrix.  The  first- 
order  step  response  is  easily  evaluated  using  the  pole-residue  informa¬ 
tion  of  the  linearized  circuit.  The  second-order  response  is  evaluated 
by  performing  the  association  of  poles  and  residues  in  a  prescribed 


manner.  The  third-order  response  is  computed  similarly.  No  numerical 
integration  is  used  in  these  computations. 

Since  the  hybrid  analysis  forms  the  basis  for  forming  the  open  cir¬ 
cuit  impedance  matrix/  the  following  linear  elements  are  allowed  by  the 
* 

program  :  resistors,  capacitors,  inductors,  voltage  or  current  sources, 
and  all  four  types  of  controlled  sources.  The  nonlinear  elements  are 
assumed  to  be  voltage  controlled,  with  the  following  polynomial  descrip¬ 
tions: 

i  =  a.fCv  3  +  a-fCv2}  +  a,fCv3]  (5-32) 

pip  c  p  j  p 

V  *  a10vq  *  a0l“r  *  a20vq  +  a02vr  * 


a.-v  v  +  a,nv 
11  q  r  30  q 


(5-33) 


where  ifl  and  vn  are  currents  and  voltages  across  branch  n,  f  is  a  linear 

^  t 

operator  of  the  type  -rr,  /,  or  constant,  and  a.,  are  constants.  It 
should  be  noted  that  eqn.  (5-33)  models  a  3-port  device. 

In  the  present  version,  PRANC  imposes  the  following  restrictions  on 
the  circuit  parameters:  maximum  number  of  elements  (both  linear  and 
nonlinear)  =  50;  maximum  number  of  nonlinear  elements  *  10;  maximum  num¬ 
ber  of  dependent  nonlinear  elements  (eqn.  5-33)  =  5;  maximum  number  of 
reactive  elements  s  20;  maximum  number  of  independent  nodes  s  30;  number 
of  input  frequencies  =  5.  These  restrictions  can  be  relaxed  if  desired. 
The  modular  structure  and  algorithms  of  PRANC  makes  it  possible  to  ex¬ 
tended  the  order  of  analysis  in  a  straightforward  manner.  The  limit  on 

*A  direct  nodal  analysis  would  only  allow  for  voltage  controlled 
current  source. 
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the  highest  order  will  eventually  be  dictated  by  the  stdrage  restric¬ 
tions  of  the  computer. 

The  validity  of  the  results  obtained  from  using  PRANC  has  been  ver¬ 
ified  through  hand-worked  examples  and  comparing  with  the  results  ob¬ 
tained  from  NCAP  C443.  In  the  next  chapter  we  will  present  two  examples 
showing  the  results  obtained  from  the  use  of  PRANC. 
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CHAPTER  6 


NONLINEAR  CIRCUIT  ANALYSIS  EXAMPLES 

6-1 .  Introduction 

In  this  section,  results  from  the  analysis  of  two  nonlinear  cir¬ 
cuits  using  PRANC  will  be  discussed.  For  the  sake  of  brevity,  the  cir¬ 
cuits  are  kept  simple.  A  brief  comparison  of  the  execution  times  for 
PRANC  and  NCAP  C443  is  also  included  in  this  section. 

Previous  works  C7, 10,363  have  used  the  Volterra  series  method  to 
study  the  distortion  phenomenon  in  transistor  amplifiers.  These  authors 
show  a  good  agreement  between  the  predicted  and  the  measured  values  of 
the  response  for  small  signal  operation.  Thus,  the  Volterra  series  is 
established  as  a  viable  approach  for  analyzing  nonlinear  circuits.  In 
section  6-2,  a  single  stage  untuned  amplifier  is  analyzed  using  PRANC. 

Time-domain  analysis  using  the  Volterra  series  method  has  not  re¬ 
ceived  much  attention  in  the  literature.  As  such  no  practical  circuits 
have  been  analyzed  using  this  technique  to  study  the  validity  of  the 
results  obtained.  In  section  6-3,  the  step  response  of  a  simple  non¬ 
linear  circuit  is  obtained.  These  results  are  then  compared  with  those 
obtained  from  a  numerical  integration  method. 

Finally,  section  6-4  compares  the  execution  times  for  the  specific 
circuit  examples  analyzed  using  PRANC  and  NCAP. 
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6-2.  Spectrum  and  Distortion  Analysis  Example. 


Consider  the  single  stage  untuned  bipolar  transistor  amplifier  cir¬ 
cuit  shown  in  Fig.  6-1.  This  circuit  configuration  has  been  taken  from 
Cl 03 .  The  results  obtained  from  using  PRANC  are  in  excellent  agreement 
with  those  obtained  from  using  NCAP,  which,  as  reported  in  [103,  are  in 
good  agreement  with  the  measured  results. 

The  equivalent  circuit  model  C73  used  for  the  bipolar  transistor  is 
shown  in  Fig.  6-2.  It  contains  three  nonlinear  elements:  a  nonlinear 
resistor,  a  nonlinear  capacitor,  and  a  dependent  nonlinearity.  The 
operating  point  of  the  transistor  is:  VC£  =  10V,  Ic  =  10mA.  The  per¬ 
tinent  parameters*  of  the  transistor  are  given  in  Table  6-1. 

Table  6-1.  Transistor  Parameters  for  the  Circuit  of  Fig.  6-1. 


Parameter 

Value 

Parameter 

Value 

11.4mA 

n 

4.6 

XC 

10.0mA 

a 

0.125 

TC 

,Cmax 

150.  mA 

8.2 

> 

10.1O 

2.30 

yFEmax 

VCB 

uCB0 

9.27  V 

re 

635kQ 

140.  V 

c  I 

1040pF 

0.348 

‘ci 

c3 

ll.lpF 

k 

25pF 

1.5pF 

The  nonlinear  elements,  using  the  Taylor  series  expansion  of  the 
nonlinear  functions,  are  represented  by  the  following  polynomials: 


♦These  parameters  depend  on  the  transistor  type 


(6-1) 


iR  =  K<V  =  91Vr  +  92VR  +  S3VR 


d  v  dvj;  dv^ 

'c  *  P(V  stcidr*  cc2  3T  ♦  CC3  3T 
S  =  S<V\b>  =  910UR  *  90l”cb  *  Wr  *  902vcb 
*  91lVcb  *  930VR  *  903Vcb  *  92lVcb  *  912VRVcb 


(6-2) 


(6-3) 


Using  the  analytical  formulae  for  g.,  C  .,  and  g..  derived  in  C7]  and 

1  Cl  1] 

the  transistor  parameters  given  in  Table  6-1,  the  numerical  values  of 
the  coefficients  of  the  polynomials  appearing  in  eqns.  (6-1)  through 
(6-3)  are  obtained.  These  values  are  given  in  Table  6-2. 

The  values  in  Table  6-2  along  with  the  linear  element  values,  the 
circuit  topology,  and  the  source  information  are  used  as  the  input  in¬ 
formation  for  PRANC.  The  effects  of  the  frequency  and  the  input  ampli¬ 
tude  on  the  output  response  are  discussed  below. 

6-2.1.  Effect  of  Frequency;  The  effect  of  frequency  on  the  system 
response  is  directly  related  to  the  effect  of  frequency  on  the  transfer 
functions.  In  the  following  paragraph  we  elaborate  on  this  concept. 

The  relationship  of  the  first-order  steady  state  response  and  the 
first-order  (linear)  transfer  function  is  well  understood.  For  the  cir¬ 
cuit  of  Fig.  6-1,  the  first  order  transfer  function  magnitude  charac¬ 
teristics  is  shown  in  Fig.  6-3.  It  should  be  apparent  that  the  amplif¬ 
ier  has  a  fairly  broadband  first  order  response. 

The  magnitude  characteristics  of  the  second-order  transfer  function 
as  a  function  of  frequency  can  be  plotted  similarly.  However,  in  this 
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Table  6-2.  Coefficients  of  the  polynomial  type  nonlinearities  in  Fig. 
6-2. 


Nonlinear 

Element 

Coefficient 

Values 

Emitter  Resistive 
Nonlinearity,  K(vR) 

g. =0.4348  - 
g^=1 .054x1 0^ 

g2=8.291 

Collector  Capacitance 
Nonlinearity,  P(vc> 

Cf1 11.1x10"^. 
C^j=4.784x10~1*> 

Cc2=-1.996x10'13 

Collector  dependent 
Nonlinearity,  G(vR,vcb) 

g-0=0.3856 

S20=7-419  -7 

g„*7.215x10.rn 

gi!  =3. 397x10  'u 
g!j2=1  -401x1 0~r 

3qi=1 *871x1 0~g 
g0P=3. 633x10  v 
g^=94.66 
g^=1 .388x10  3 
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case,  two  frequencies  are  involved,  and  the  significance  of  the  magni¬ 
tude  plot  must  be  clearly  understood.  An  example  of  this  plot  is  shown 
in  Fig.  6-4,  which  depicts  the  magnitude  characteristic  of  the  second 
harmonic  frequency  response  as  the  input  tone  is  varied  from  lOOKHz  to 
1 00MHz .  Similarly,  if  the  magnitude  characteristics  of  an  intermodula¬ 
tion  frequency,  wIM,  were  desired,  a  similar  plot  of  .,jw^) | ,  such 

that  u>.  +  w.  =  in  the  band  of  interest,  would  provide  the  informa- 
i  k  IM 

tion.  Thus,  by  letting  uk  =  |wIM  “  and  sweeping  across  the  band 
of  interest,  the  information  on  how  the  magnitude  of  an  intermodulation 
product  will  vary  can  be  obtained. 

The  magnitude  characteristic  of  the  third-order  transfer  function 
has  a  similar  significance.  In  this  case,  a  combination  of  three 
discrete  frequencies  is  involved. 

6-2.2.  Effect  of  Input  Amplitude.  The  effect  of  the  input  amplitude 
on  the  distortion  products  for  the  amplifier  circuit  is  shown  in  Fig. 
6-5.  A  single  tone  at  2.5MHz  is  used  as  the  input  to  the  amplifier  and 
the  second  and  third  harmonic  outputs  are  computed.  On  the  double  loga¬ 
rithmic  paper,  we  note  that  the  distortion  plots  are  straight  lines, 
with  m3  >  m^  >  where  m^,  n^,  and  m^  are  the  slopes  of  the  first, 
second,  and  third  harmonic  distortion,  respectively.  It  is  noted  that 
the  validity  of  the  small-signal  model  extends  up  to  approximately  25mV 
input  amplitude,  after  which  the  Volterra  series  prediction  is  more  pes¬ 
simistic  than  the  measured  results  for  the  amplifier  C1Q3.  Also,  the 
harmonic  distortion  indices  for  a  given  input  amplitude  can  be  deter- 

*Since  FRANC  computes  the  entire  output  spectrum,  it  automatically 
takes  the  negative  frequencies  into  account.  Thus  the  user  does  not 
have  to  specify,  for  example,  a  -  2.5MHz  and  3MHz  frequency  to  get 
the  intermodulation  product  at  0.5MHz. 
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Output  (Vo)ts) 


Input  (Volts) 

Figure  6-5.  Harmonic  Distortion  Variation  With  Input  Amplitude 


mined  directly  from  Fig.  6-5. 

6-3.  Time-Domain  Analysis  Example.  Consider  the  simple  resistive  non¬ 
linearity  circuit  of  Fig.  6-6.  The  v-i  characteristics  for  the  non¬ 
linear  element  is  i=10v^.  In  this  section  we  will  determine  the  step 
response  of  the  circuit  using  Volterra  series  method  and  then  compare 
the  results  obtained  from  an  exact  numerical  integration  method. 


By  inspection,  the  first-order  transfer  for  the  circuit  is: 


W  *  '  r  *  we- » '  -r-  “  =  i'  Gi  *  iq’  <6-‘> 

For  an  input  Au(t),  the  first  order  response,  y^Ct),  is  given  by 

y1 <t>  =  a(1  -  e”pt)  ;  a  =  ~-  (6-5) 


Since  there  is  no  square  term  in  the  nonlinearity,  the  second-order 
transfer  function  is  zero.  The  third-order  transfer  function  is  found 


-  3exp(-2pt)  +  0.5exp(-3pt)3  (6-9) 

The  response  y(t)  can  be  calculated  using  eqn.  (6-9).  The  step  response 
of  the  circuit  of  Fig.  6-8  is  not  as  interesting  as  the  error  incurred 
in  using  Volterra  series.  A  summary  of  the  error  vs.  the  l^  norm  of  the 
linear  kernal  function  is  given  in  Table  6-3.  The  norm  is  varied  by 
changing  the  values  of  in  the  circuit.  This  alters  the  value  of  p, 
but  not  of  r,  in  eqn.  (6-4).  An  estimate  of  the  truncation  error  as 
predicted  by  Theorem  2-2  is  also  given  in  Table  6-3.  Clearly,  there  is 
a  good  agreement  between  the  calculated  and  the  predicted  values  of  this 
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error. 


Table  6-3.  Summary  of 

Truncation  Error 

l.j  Norm  of 

Predicted  Error* 

Calculated 

Error 

h(t) 

0(30Z65X5> 

Maximum 

Steady 

State 

0.1333 

io"2 

io'2 

io'2 

0.0909 

10~3 

io-3 

10‘3 

9.9x10-3 

10‘9 

10~5 

io”10 

★"Order  of  " 

6-4.  PRANC  vs, 

■  NCAP 

Both  PRANC 

and  NCAP  C443  are  based 

on  the  Volterra 

series  method 

for  analyzing  nonlinear  circuits.  The  basic  algorithms  used  by  the  two 

programs  are,  however,  different.  PRANC  uses  a  semi-symbolic  analysis 

procedure,  where  as  NCAP  uses  a  purely  numerical  approach.  Furthermore, 

the  features  offered  by  the  two  programs  are  different.  NCAP  allows  a 

free-format  input  and  computes  the  transfer  functions  and  node  voltages 

at  2n-1  discrete  frequency  points,  where  n  is  the  order  of  the 

analysis.*  PRANC,  on  the  other  hand,  requires  a  fixed  format  input,  but 

computes  the  complete  output  spectrum.  Thus,  the  number  of  output  fre- 

*For  NCAP,  this  happens  to  be  the  number  of  input  frequencies  also, 
since  the  user  must  specify  n  input  tones  in  order  to  obtain  an  n- 
order  analysis. 
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quency  points  for  PRANC  is  much  higher  than  that  for  NCAP.  Because  of 
the  algorithm  and  the  feature  differences,  the  execution  times  required 
by  the  two  programs  for  analyzing  similar  circuits  are  different.  In 
the  following  paragraph  we  compare  the  approximate  execution  times  for 
the  two  programs  on  the  COC  6500  computer  at  Purdue  University. 

The  bipolar  transistor  amplifier  of  Figure  6-1  was  analyzed  on  both 
PRANC  and  NCAP.  For  a  two-tone  input  and  up  to  a  third-order  analysis, 
PRANC  required  1 .74  sec  to  computed  the  complete  output  spectrum  (in¬ 
volving  approximately  32  frequency  points).  For  a  5-tone  input,  it  re¬ 
quired  6  seconds  to  do  the  same  job  (approximately  150  points).  By  con¬ 
trast,  NCAP  required  7.5  sec  to  compute  the  response  at  7  discrete  fre¬ 
quencies  for  a  three  tone  input.  Other  examples  exercised  on  PRANC  and 
NCAP  show  a  similar  advantage  in  execution  times:  PRANC,  despite  pro¬ 
viding  more  information  to  the  user,  requires  less  execution  time  then 
NCAP.  A  more  detailed  comparison  of  the  two  programs  is  contained  else¬ 
where  C413. 
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CHAPTER  ? 


CONCLUDING  REMARKS 


7-1 .  Summary 

The  contents  of  this  thesis  has  dealt  with  the  various  aspects  of 
Volterra  series  analysis  of  lumped  systems,  with  polynomial  type  non¬ 
linear  elements,  with  special  emphasis  on  nonlinear  circuits.  The  fun¬ 
damental  intent  behind  this  work  is  to  show  that  the  Volterra  series 
method  could  be  used  in  a  variety  of  engineering  calculations  -  particu¬ 
larly  computer-aided  analysis  of  mildly  nonlinear  systems. 

By  investigating  a  concrete,  ye..  general,  set  of  differential  equa¬ 
tion  that  characterize  lumped  nonlinear  systems,  a  unified  theory  has 
been  developed  here  which  can  be  used  to  study  the  behavior  of  such  sys¬ 
tems. 

By  analyzing  how  the  series  solution  converges,  simple  and  easy  to 
use  recursive  relationships  have  been  derived.  These  relationships  can 
be  used  for:  1)  estimating  the  error  incurred  in  truncating  the  series; 
and  2)  determining  the  bound  on  the  input  function  for  which  the  series 
converges.  The  higher-order  kernel  functions  are  not  required  a  priori 
in  these  relationships;  only  the  linearized  system  need  be  analyzed  to 
determine  how,  if  at  all,  the  series  converges. 

The  determination  of  the  transfer  function  by  the  direct  applica¬ 
tion  of  the  multi-dimensional  transforms  h&s  been  developed  here.  This 
approach  provides  a  more  direct  and  rigorous  way  of  characterizing  the 
system  using  Volterra  series. 

The  entire  investigation  culminated  in  the  development  of  a  comput¬ 
er  program  which,  through  the  use  of  semi-symbolic  analysis,  provides  an 
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efficient  tool  for  spectrum  and  distortion  analysis  of  mildly  nonlinear 
circuits.  The  time-domain  analysis  of  nonlinear  circuits  has  also  been 
investigated. 

7-2.  Further  Research 

There  are  several  problems  that  were  unraveled  during  the  course  of 
this  investigation.  We  mention  a  few  of  them  here. 

In  dealing  with  the  convergence  of  a  class  of  lumped  nonlinear  sys¬ 
tems,  we  showed  that  the  truncation  error  involved  three  quantities:  1) 
the  l ^  norm  of  the  linearized  system,  2)  the  bound  on  the  system  input, 
and  3)  the  coefficients  of  the  polynomial  describing  the  nonlinear  ele¬ 
ments.  An  important  consequence  of  this  result  leads  to  the  following 
question:  How  do  we  speed  up  the  convergence  of  the  series  solution  and 
yet  obtain  meaningful  results?  One  factor  which  has  a  direct  bearing  on 
the  answer  to  this  question  is  the  criterion  and  method  used  for  approx¬ 
imating  the  nonlinear  function  by  polynomials.  By  changing  the  coeffi¬ 
cient  values,  the  convergence  would  also  change.  These  coefficients 
depend  on  the  norm  of  the  error  function  we  choose  to  minimize  when  ap¬ 
proximating  the  nonlinear  function  by  a  polynomial.  Thus,  the  need  for 
establishing  a  criterion  (choice  of  a  norm)  to  arrive  at  a  rapidly  con¬ 
verging  series  is  important. 

Another  problem  is  related  to  the  interpretation  of  the  nonlinear 
transfer  functions.  We  have  shown  that  these  can  be  obtained  directly 
from  the  application  of  multidimensional  transforms.  In  the  case  of 
linear  systems,  the  poles  and  zeros  have  a  direct  relation  to  the  system 
response  and  stability.  Since  the  Volterra  series  solution  regards  the 

linear  system  as  a  limiting  case  of  the  nonlinear  system,  an  important 
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question  to  ask  is:  Is  there  a  relation  between  the  nonlinear  transfer 
functions  and  the  system  stability?  How  does  one  interpret,  for  exam¬ 
ple,  the  significance  of  in  the  (s^+s^+P^)  term  in  the  second-order 
transfer  function?  If  these  questions  can  be  adequately  answered  in 
terms  of  the  circuit  parameters,  a  significant  insight  into  the  behavior 
of  nonlinear  systems  would  be  gained.  As  pointed  out  earlier,  working 
in  the  transform  domain  is  much  simpler  and  algorithmic  as  compared  to 
working  in  the  time  domain. 

The  synthesis  problem  using  Volterra  series  has  not  received  much 
attention.  It  is  an  important  problem,  one  which  can  be  used  in  a 
variety  of  applications.  One  of  them,  for  example,  could  be  to  improve 
the  received  signal-to-noise  ratio  of  a  signal  subjected  to  nonlinear 
processing. 
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APPENDIX 

A.  For  Chapter  2 


A-1 .  Derivation  of  Eqn.  2-23; 

According  to  eqn.  (2-4); 

y2<t)  *  h2<T^T2^x^t“T‘j)x{t”T2>dTidT2  (A-1) 

From  eqn.  (2-19),  we  have 

y2(t)  a  -«a2  US  h(t-t)R (T-Ti )R<t-t2) - 


xCyxC^Jd^d^d* 

We  make  the  following  change  of  variables: 


(A-2) 


u 


=  t  -  T 
=  t  -  T 


1 


S>  T 
=>  T, 


u2  *  '  -  ’2  “> 


T-  =  t  -  T_ 
2  2 


Substituting  eqn.  (A-3)  in  eqn.  (A-2),  we  get 


(A-3> 


y2(t)  a  -ea2  ss:s  h(u)  R(u^-u)  fIXUg-u) 


y 2 Ct )  =  -ea^  h(x)  fi(Ti  "  T>  R(t2  "  t)* 


x(t  -  T  )  x(t  “  x  )dTdt.dT_ 
1  2  12 


(A-5) 


Comparing  (A-5)  with  (A-1),  we  get 


*n2('X‘\'X2*  ~  ta  iSo  h(x)  R(t^-t)  R(T2“T)dT 


(A-6) 


A-2.  Proof  of  Theorem  2^2} 
Given 


L(p)  y (t)  +  £  a  yn<t)  =  x(t) 
_■»  n 
n=2 


Applying  the  inverse  operator  we  get 


y<t)  +  /o  h<t  "  T) 


£  a  yn(T)  dT  = 

|n-2  n 


JQ  h(t  -  T)x(t)dT  (A-7) 


Assuming  a  series  solution  for  y(t),  such  that 


y(t)  =  y^  (t)  +  y2(t)  +  y3(t>  +  ••• 


(A-8) 


and  substituting  in  eqn.  (A-7),  we  get 


£  y . (t)  +  f*  h(t 
i=1  J0 


[„4  *"  'I,  ”<‘,>,n] dT  ■  •£ 


h(t  -  t)  x(t)dx 


(A-9) 


Solving  for  y,j(t),  y2(t),  ...  successively  In  eqn.  (A-9),  we  get 
the  following  set  of  equations: 


y.j<t)  -  hCt  *  x)  x(t)  dt 

y2<t)  +  a £  h(t  -  t)  y^(t)  dr  =  0 

y3(t)  +  a2  h<t  -  t)  Ja2)'1  (Oyjjtt)  +  ^y^^J 

y.(x)  +  f  h(t  -  t)  i  (t)  dx  =  0 

4  *0  n 

where  i  (t)  is  shown  CIO-123  to  be: 
n 

n 


i  <t>  =  £  aky  . 
n  lcs2  ^ 

L-m+1 


y,  m  =  Z  Y.Y.  .  - 

Then,  for  G  =  |h(t)|dt  and 


X  =  max  [x(t>3 

0  <  t  <  • 


n  >  2 


dt 


0 


(A-10) 


(A-11) 


we  get 


Y^  *  |y(t) |  <  GX 

y2  1  la2l<i3x2 

Y  <  G|I  | 
n  -  1  n' 

where  Yn  is  given  by  eqn.  (2-38)  after  we  use  eqn.  (A-11).  The  inequal¬ 
ity  (2-37)  then  follows  fro«  eqn.  (A-8). 
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B.  For  Chapter  3 


B-1 .  Derivation  of  *C£-(t) 3 : 

Using  the  functional  property  of  the  6-function,  we  can  write 

y2(t)  *  /Lf  y<T1)y1<T2)S<t  -  T.,)6(t  -  t2>  d^d^  (B-1) 

Introducing  dummy  variables  t^  and  t2  in  eqn.  (B-1),  we  get 

y(trt2)  =  y2(t)  =  f_J  y(T1)y(T2)fi(t1  -  T1>6(t2  -  t?)  d*^  (b-2) 

Taking  a  two-dimensional  transform  of  eqn.  (B-2),  as  per  eqns. 
(3-11)  and  (3-12),  we  get 

y(s.,,s2)  =  ‘[^VV]  =  y(si>y(s2)  <B-3) 

The  derivation  for  ljyn(t)j  is  done  similarly  by  introducing  n-delta 
functions  first,  etc. 

B-2.  Derivation  of  yZ(t)J 

Using  the  chain  rule,  we  get 

~y2(t)  =  2y(t)  y(t)  =  2y(t>  y'(t)  (B-4) 

Using  eqn.  (B-2),  we  re-write  eqn.  (B-4)  as: 

2y(t>  y'M  [2  X<«,>  x'<'2)  ««,  -  «,)  «<tj  -  t2>  d^dTj 

y<T2)y'<t1)  «(t1  -  )  6(t2  -  t2>  dt^Tg  (B-5) 

Taking  a  two  dimensional  transform  of  eqn.  (B-5),  we  get 
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(B-6) 


Yls^Sg) 


=  l[: 


2s2Y(s1)Y(s2) 


2s1Y(s1)Y(s2) 


=  (s1  +  s2)Y(s1>Y(s2) 


The  derivation  of  i 


is  done  similarly. 


